2.1. Governing equations

Figure 1.

Schematic illustration of the two-dimensional flow domain. The regular porous media is represented by the grey obstacles. The size of a basic unit is $D^*=d^*+l^*$ .

We consider a 2-D RB system filled with regular square obstacles, as shown in figure 1. Constant species concentrations are kept on the two horizontal plates with a separation of $H^*$ . Hereafter, the asterisk $*$ represents the dimensional forms of the variables; $c^*$ is the relative species concentration, which is equal to zero at the bottom plate and $\varDelta _c^*$ at the top plate. The fluid density is then determined by $c^*$ as $\rho ^*=\rho _0^*(1+\beta c^*)$ , in which $\rho _0^*$ is the density of a reference state and $\beta$ is the expansion coefficient of the concentration. The periodic boundary condition is used in the horizontal direction of the flow domain. The no-slip boundary condition is applied to all solid surfaces, including the two plates and the fluid–obstacle interfaces. In addition, the obstacles are assumed to be non-penetrated by the species. Within the Oberbeck–Boussinesq approximation, the incompressible flow in the pores is governed by

(2.1a) \begin{eqnarray} \frac {\partial \boldsymbol {u}^*}{\partial t^*} + \boldsymbol {u}^*\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u}^* &=& -\boldsymbol {\nabla } p^* + \nu \boldsymbol {\nabla }^2\boldsymbol {u}^* - g \beta c^*\boldsymbol {e}_z, \end{eqnarray}

(2.1b) \begin{eqnarray} \frac {\partial c^*}{\partial t^*}+ \boldsymbol {u}^*\boldsymbol {\cdot }\boldsymbol {\nabla } c^* &=& \kappa \nabla ^2 c^*, \end{eqnarray}

(2.1c) \begin{eqnarray} \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u}^* &=& 0. \end{eqnarray}

Here, $\boldsymbol {u}^*$ is the velocity vector, $p^*$ is pressure, $\nu$ is viscosity, $g$ is the gravitational acceleration, $\boldsymbol {e}_z$ is the vertical unit vector and $\kappa$ is the mass diffusivity of the species. We use the domain height $H^*$ , the free-fall velocity $\sqrt {gH^*\beta \varDelta _c^*}$ and the concentration difference $\varDelta _c^*$ to non-dimensionalise the governing equations (2.1). Dimensionless control parameters include the aspect ratio, the Schmidt number and the Rayleigh number, which are defined respectively as

(2.2) \begin{equation} \Gamma = \frac {L^*}{H^*}, \quad {\ {Sc}} = \frac {\nu }{\kappa }, \quad {\ {Ra}}=\frac {g\beta \varDelta _c^* {H^*}^3}{\kappa \nu }. \end{equation}

Here, $L^*$ is the domain width. Then, the dimensionless governing equations read

(2.3a) \begin{eqnarray} \frac {\partial \boldsymbol {u}}{\partial t} + \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }\boldsymbol {u} &=& -\boldsymbol {\nabla } p + \sqrt {\frac {{\ {Sc}}}{{\ {Ra}}}}\boldsymbol {\nabla }^2\boldsymbol {u} - c\boldsymbol {e}_z, \end{eqnarray}

(2.3b) \begin{eqnarray} \frac {\partial c}{\partial t}+ \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla } c &=& \frac {1}{\sqrt {{\ {Sc}}{\ {Ra}}}} \nabla ^2 c, \end{eqnarray}

(2.3c) \begin{eqnarray} \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {u} &=& 0, \end{eqnarray}

and the boundary conditions at two plates read

(2.4a) \begin{eqnarray} \boldsymbol {u}&=&0, \quad c=0,\qquad \text{at}\ z=0; \end{eqnarray}

(2.4b) \begin{eqnarray} \boldsymbol {u}&=&0, \quad c=1,\qquad \text{at}\ z=1. \end{eqnarray}

The key response parameters include the Sherwood number and the Reynolds number, which are defined respectively as

(2.5) \begin{equation} {\ {Sh}} = \frac {\kappa \left\langle \frac {\partial c^*}{\partial z^*}\right\rangle _{z^*=0} }{\kappa \varDelta ^*_c {H^*}^{-1}}=\left\langle \frac {\partial c}{\partial z}\right\rangle _{z=0} , \quad \ {Re} = \frac {\boldsymbol {u}^*_{rms} H^*}{\nu }=\sqrt {\frac {{\ {Ra}}}{{\ {Sc}}}}\boldsymbol {u}_{rms}. \end{equation}

Hereafter, the bracket $\langle \cdot \rangle$ denotes the average over a specific horizontal plane. The subscript ‘rms’ stands for the root-mean-square value over the entire domain except for obstacles. Note that the Sherwood number $\ {Sh}$ represents the ratio of the total mass transfer rate to the rate of diffusive mass transport. In the statistically steady state, $\ {Sh}$ should be the same at the two plates.

2.2. Porous media

The porous medium structures set in this study are all regularly arranged squares, as shown in figure 1. The side length of each square is $d^* (=dH^*)$ and the spacing between them is $l^* (=lH^*)$ . In traditional macroscopic simulations of convection in porous media, the REV which contains several pores is generally defined as the basic research unit (Nield & Bejan Reference Nield and Bejan2017). In our pore-scale simulations, the basic unit consists of an obstacle and its surrounding void space, whose side length is $D^*=d^*+l^*$ , as shown in figure 1. Under this definition, the whole domain is filled with the basic units. The additional control parameters related to the porous medium structure include the porosity, the Darcy number and the Rayleigh–Darcy number, which are defined as follows:

(2.6) \begin{equation} \phi =1-\left(\frac {d^*}{D^*}\right)^2,\quad {\ {Da}}=\frac {K^*}{{H^*}^2}, \quad {\ {Ra}}_D=\frac {g\beta \varDelta _c^* H^* K^*}{\kappa \nu }={\ {Ra}} {\ {Da}}. \end{equation}

Here, $K^*$ is the permeability, which can be calculated by Darcy’s law (Nield & Bejan Reference Nield and Bejan2017):

(2.7) \begin{equation} K^*=-\nu \frac {u_m^*}{\nabla p*}. \end{equation}

Here, $u_m^*$ is the mean velocity of the flow induced by the pressure gradient $\nabla p^*$ . Note that $\nu$ represents the kinematic viscosity and $p^*$ includes the density. Table 1 summarises the three sets of porous medium configurations that we have defined, with porosity of 0.64, 0.36 and 0.15, respectively. In these configurations, we keep the size of the obstacles constant and vary their number to change the porosity. Accordingly, the spacing between the obstacles and the size of the basic unit also change.

Table 1.

Details for the porous media. Columns from left to right are: the porosity, the aspect ratio, the number of obstacles in the horizontal direction and the vertical direction, the size and the spacing of the obstacles, the size of the basic unit, the two Darcy numbers calculated by the Darcy’s law, and Kozeny’s equation with $\eta =125$ .

To calculate the permeability, we conducted additional numerical simulations for each porosity. In these simulations, the effect of the concentration field is neglected and a uniform horizontal pressure gradient is applied throughout the domain. The horizontal mean velocity is then calculated after the flow stabilises. During the averaging process, we exclude regions within five basic unit distances near the upper and lower boundaries to minimise boundary effects. Note that the current porous media structure consists of obstacles arranged in a regular pattern. Therefore, the permeability in both the horizontal and vertical directions should be the same, and it is sufficient to calculate only one of them. Additionally, we ensure that the pore-scale Reynolds number, i.e.

(2.8) \begin{equation} \ {Re}_p = \frac {|\boldsymbol {u}^*| l^*}{\nu }=l\sqrt {\frac {{\ {Ra}}}{{\ {Sc}}}}|\boldsymbol {u}|, \end{equation}

for these cases is much less than 1, so the flow remains in the Darcy regime (Hewitt Reference Hewitt2020).

Figure 2.

Flow driven solely by pressure gradient. (a) Horizontal mean velocity $u_m$ versus the horizontal pressure gradient $-\nabla _x p$ . (b) Darcy number $Da$ versus the porosity $\phi$ . In panel (b), the values represented by the symbols are calculated from the slope of the corresponding fitted line in panel (a), and the dashed line is derived from Kozeny’s equation (2.10) with $\eta =125$ .

Figure 2(a) illustrates the relationship between the horizontal mean velocity $u_m$ and pressure gradient $\nabla _x p$ for different porosities. It is evident that as $-\nabla _x p$ increases, $u_m$ increases linearly. By substituting (2.7) into (2.6), we can derive the formula for calculating the Darcy number:

(2.9) \begin{equation} {\ {Da}}=-\sqrt {\frac {{\ {Sc}}}{{\ {Ra}}}}\frac {u_m}{\nabla _x p}. \end{equation}

In the above cases, we set ${\ {Sc}}=400$ and ${\ {Ra}}=10^6$ . By performing a linear fit on the data points in figure 2(a), we obtain the slopes, which can then be used to calculate the Darcy number using (2.9). The results are presented in figure 2(b). It can be observed that the Darcy number increases nonlinearly with porosity. It should be noted that the permeability can also be estimated by Kozeny’s equation (Nield & Bejan Reference Nield and Bejan2017):

(2.10) \begin{equation} K^*=\frac {{d^*}^2\phi ^3}{\eta (1-\phi )^2}, \end{equation}

in which $\eta$ is an empirical coefficient. The classical Carman–Kozeny equation provides $\eta =180$ , which applies to media composed of regularly packed spherical particles. The Blake–Kozeny equation gives $\eta =150$ , suitable for more complex configurations (Dullien Reference Dullien2012). In the current study, $\eta =125$ aligns well with the simulation data, as indicated by the dashed line in figure 2(b). This value is nearly identical to $\eta =126$ obtained in the simulations of Gasow et al. (Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020), where their porous medium configuration also consisted of regularly arranged squares. The specific values of the Darcy number are also summarised in table 1. It should be noted that in actual CO $_2$ sequestration, the Darcy number is very small, typically reaching the order of $10^{-13}$ (Jin & Kuznetsov Reference Jin and Kuznetsov2024). For pore-scale simulations, achieving such small pore sizes would require very high resolution, which is computationally infeasible with current capabilities. In the series of works by Gasow et al. (Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020, Reference Gasow, Kuznetsov, Avila and Jin2021, Reference Gasow, Kuznetsov and Jin2022), the lowest Darcy number achieved was of the order of $10^{-8}$ , the same as in this study. Nevertheless, the flow behaviours observed at larger Darcy numbers show a certain consistency, allowing for reasonable extrapolation to cases with smaller Darcy numbers.

In the numerical simulations of this study, the boundary conditions of the obstacles are implemented using the direct-forcing immersed boundary method (IBM) (Uhlmann Reference Uhlmann2005). Specifically, an additional body force $\boldsymbol {f}$ is introduced in the momentum equation (2.3a ) and the convection–diffusion equation (2.3b ) to ensure zero velocity and zero normal gradient of concentration at the obstacle boundaries. Particularly, at the corners of the square, the sum of the fluxes of the concentration field along the $x$ and $z$ directions is set to zero. Under such settings, the flow state inside the obstacles does not affect the external flow field. This method has been widely used in various complex geometries and deformable interfaces (Spandan et al. Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017; Vanella & Balaras Reference Vanella and Balaras2020), and was also employed in our previous simulations involving circular obstacles (Liu et al. Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020). For the current square obstacles, their boundaries are directly defined on the Eulerian grid, thus avoiding additional errors due to interpolation. To ensure accuracy, we extend two grid layers inward from the obstacle boundaries, enforcing both no-slip and no-penetration boundary conditions. Additionally, we have conducted some a posteriori validations, with details provided in the next section.

In the current simulations, the obstacles are arranged in a regular pattern, which differs from the real irregular porous structures found in geological formations. However, it is extremely challenging to set up randomly distributed obstacles in pore-scale simulations, especially in cases with very low porosity. The difficulty mainly lies in the narrow spacing between obstacles caused by randomness and the complex structures formed by the combination of multiple obstacles. These require highly refined grids for resolution, which poses a significant challenge given the current computational capabilities. Therefore, we continue to focus on cases with regularly arranged obstacles. We believe that the conclusions drawn from this basic configuration can provide us with a deeper understanding of the physical mechanisms of convection in porous media. Additionally, we conducted a separate simulation with randomly distributed obstacles for a porosity of 0.64, and the results are summarised in Appendix A. It can be observed that the qualitative trend of mass transport is consistent between the regular and random distributions of obstacles, although there are some minor quantitative differences. A systematic study on this aspect is still needed in future research.

2.3. Numerical settings

We conduct direct numerical simulation (DNS) by solving the governing equations (2.3), together with the additional body force $\boldsymbol {f}$ mentioned above. The code used in this study is an improved version based on our in-house code, which has been widely used in research on double-diffusive convection (Yang et al. Reference Yang, van der Poel, Ostilla-Mónico, Sun, Verzicco, Grossmann and Lohse2015, Reference Yang, Verzicco, Lohse and Caulfield2022; Li & Yang Reference Li and Yang2023). Basically, the finite difference method and the fractional time-step method are employed (Ostilla-Mónico et al. Reference Ostilla-Mónico, Yang, van der Poel, Lohse and Verzicco2015). One advantage of this code is the use of a dual-resolution technique, where a more refined grid is employed for the scalar field with a high Schmidt number $\ {Sc}$ . In our previous studies, $\ {Sc}$ reached as high as 1000 (Yang et al. Reference Yang, Verzicco, Lohse and Caulfield2022). In the current study, we fix ${\ {Sc}}=400$ , a value close to the practical conditions of CO $_2$ sequestration (Huppert & Neufeld Reference Huppert and Neufeld2014). For the three types of porous medium configurations shown in table 1, we have set a wide range of Rayleigh number $\ {Ra}$ , with the corresponding Rayleigh–Darcy number ${\ {Ra}}_D$ ranging from $10^2$ to $10^6$ . In the classical Darcy model, ${\ {Ra}}_D \gtrsim 1300$ signifies the high Rayleigh regime (Nield & Bejan Reference Nield and Bejan2017). Therefore, our parameter range extends far beyond this regime to study the transition from the Darcy regime to the non-Darcy regime. The details of the simulations are summarised in table 2.

Table 2.

Numerical details for all the cases. Columns from left to right are: the porosity, the Rayleigh number, the Rayleigh–Darcy number, the aspect ratio, the resolutions with refined factors in the horizontal and vertical directions, the simulation time before the statistical stage, the statistical time, the statistical Sherwood number, the statistical Reynolds number, and the relative difference of the statistical Sherwood numbers calculated at the two plates.

It should be noted that the aspect ratio $\Gamma$ of the domain can also influence flow structures and mass transport, especially when the aspect ratio is very small, as the boundaries may constrain the development of convection. This phenomenon has been extensively studied in classical RB convection (van der Poel et al. Reference van der Poel, Stevens, Sugiyama and Lohse2012; Wagner & Shishkina Reference Wagner and Shishkina2013; Shishkina Reference Shishkina2021). In the current study, since our focus is not on investigating the effects of aspect ratio, we aim to set a sufficiently large $\Gamma$ within the limits of available computational resources to minimise its impact on the flow. For most cases with $\phi =0.64$ and $0.36$ , we set the aspect ratio $\Gamma =2$ ; for cases with $\phi =0.36$ and ${\ {Ra}}=10^{12}$ , as well as all cases with $\phi =0.15$ , we set $\Gamma =1$ to save computational resources. In the present configuration, the cases with the lowest $\ {Ra}$ for each porosity include at least two pairs of convection cells in the flow field. Furthermore, as we will demonstrate later, the results across cases with different aspect ratios exhibit consistency. Therefore, we consider the current aspect ratios to be acceptable.

The grid resolution meets the following three criteria: first, the gap between two obstacles must be at least five grid cells, which implies that the base grid scale is less than $0.2l$ ; second, the base grid scale is smaller than the Kolmogorov scale $\eta _K=(\nu ^3/\epsilon )^{1/4}$ (Grötzbach Reference Grötzbach1983), in which $\epsilon$ is the mean viscous dissipation rate; third, the refined grid scale is smaller than the Batchelor scale $\eta _B=(\nu \kappa ^2/\epsilon )^{1/4}$ (Silano et al. Reference Silano, Sreenivasan and Verzicco2010). For all cases, we simulated for a sufficiently long time $t_d$ to ensure that the flow reached a statistically steady state. Then, we continued the simulation for a period $t_s$ for statistical analysis. The relative difference between the time-averaged Sherwood number $\overline {{\ {Sh}}}$ in the first and second halves of this period did not exceed $1\,\%$ . Hereafter, the overline denotes time-averaged values. Additionally, we calculated the relative difference of $\overline {{\ {Sh}}}$ at the upper and lower plates, which was less than $4\,\%$ for all cases. This indicates that no additional mass flux entered or exited the obstacles.

3.1. Transition of the flow field

Figure 3.

Snapshots of the instantaneous concentration fields for the cases with (a) ${\ {Ra}}=10^7$ , (b) ${\ {Ra}}=2 \times 10^8$ , (c) ${\ {Ra}}=2 \times 10^9$ and (d) ${\ {Ra}}=10^{11}$ . The porosity is fixed at 0.64. The obstacles are denoted by the white squares.

In this section, we examine the characteristic structures of the flow field under varying $\ {Ra}$ or ${\ {Ra}}_D$ . Based on the Darcy model, convection in 2-D porous media can be broadly classified into three regimes (Hewitt et al. Reference Hewitt, Neufeld and Lister2014; Nield & Bejan Reference Nield and Bejan2017; Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020): the conducting regime ( $0 \leqslant {\ {Ra}}_D \leqslant 4\pi ^2$ ), where no flow occurs; the quasi-steady regime ( $4\pi ^2 \leqslant {\ {Ra}}_D \lesssim 1300$ ), characterised by periodically distributed large-scale convection rolls; and the high Rayleigh regime ( ${\ {Ra}}_D \gtrsim 1300$ ), where the large-scale convection rolls vanish and are replaced by irregularly distributed mega-plumes in the central region and proto-plumes near the boundaries. Figure 3 presents four typical concentration fields for the cases of $\phi =0.64$ . When ${\ {Ra}}=10^7$ and ${\ {Ra}}_D=2.6 \times 10^2$ , the flow field is composed of four pairs of convection cells, corresponding to the quasi-steady regime, as shown in figure 3(a). At this stage, the concentration distribution in the region near the upper and lower boundary plates is very uniform over at least one basic unit. As the Rayleigh number increases, large-scale structures in the flow field begin to become disordered, entering the high Rayleigh regime, as shown in figure 3(b). At this stage, large-scale convection cells disappear, replaced by smaller-scale plumes that are chaotically distributed. These plumes originate from the upper and lower boundaries, and merge into larger-scale plumes in the middle of the domain. When ${\ {Ra}}_D$ further increases to $5.2 \times 10^4$ , as shown in figure 3(c), the characteristic structures of the concentration field comprise two parts: one consists of two pairs of large-scale convection cells with a fan-like shape probably caused by mechanical dispersion (Wen et al. Reference Wen, Chang and Hesse2018), and the other consists of small-scale plumes generated at the boundaries, with sizes smaller than the gaps between obstacles. Finally, for the highest Rayleigh number in this set of cases, i.e. ${\ {Ra}}=10^{11}$ and ${\ {Ra}}_D=2.6 \times 10^6$ , the bulk area is thoroughly mixed, with clearer convection cells, as shown in figure 3(d). Meanwhile, the extremely small-scale plumes near the boundaries become highly turbulent. The flow field at this stage resembles classical RB convection, with obstacles having minimal impact.

Figure 4.

Snapshots of the instantaneous pore-scale Reynolds number $\ {Re}_p$ for the cases with (a) ${\ {Ra}}=10^7$ , (b)  ${\ {Ra}}=2 \times 10^8$ , (c) ${\ {Ra}}=2 \times 10^9$ and (d) ${\ {Ra}}=10^{11}$ . The porosity is fixed at 0.64. The obstacles are denoted by the white squares.

It is important to note that the chaotic structures in figure 3(b) differ from classical turbulence, with a very low Reynolds number of approximately 2 (see table 2). In the field of porous media convection, this type of flow is often referred to as ‘pseudo-turbulence’ (Jin & Kuznetsov Reference Jin and Kuznetsov2024). By defining the pore-scale Reynolds number $\ {Re}_p$ (2.8), we can better characterise the turbulent flows in porous media. A Darcy-type flow occurs when $\ {Re}_p\ll 1$ (Hewitt Reference Hewitt2020; De Paoli Reference De Paoli2023). Figure 4 presents the contour plots of $\ {Re}_p$ corresponding to four instantaneous fields shown in figure 3. It can be observed that in the first two cases, the magnitude of $\ {Re}_p$ is much less than 1, indicating that the flow in the pores is not turbulent and conforms to Darcy flow. However, for the flow field in figure 4(b), the larger-scale structures are random and chaotic, leading us to classify this as ‘pseudo-turbulence’. For the cases in figures 4(c) and 4(d), the maxima of $\ {Re}_p$ approach or exceed 1, clearly indicating that the classical Darcy model is no longer applicable at this stage. Furthermore, we can examine the transition process from Darcy flow to non-Darcy flow through the distribution of $\ {Re}_p$ . At ${\ {Ra}}_D=2.6 \times 10^2$ , the flow field consists of relatively stable convection cells, with alternating regions of higher and lower velocities in the bulk. When ${\ {Ra}}_D=5.2 \times 10^3$ , the chaotic distribution of $\ {Re}_p$ clearly represents the pseudo-turbulent state at this stage. At ${\ {Ra}}_D=5.2 \times 10^4$ , two pairs of large-scale convection cells form, with higher velocities primarily appearing at the boundaries of the cells, while the velocities in their centres are very low, as seen in the dark areas of figure 4(c). At this stage, some smaller-scale plumes still appear in the bulk region. For ${\ {Ra}}_D=2.6 \times 10^6$ , the bulk region is essentially dominated by only two pairs of convection cells, with very low velocities in the large central areas of these cells.

3.2. Horizontal characteristic scale

We have observed that as the Rayleigh number increases, the flow state gradually transitions from the Darcy flow to the non-Darcy flow, and the transitional state exhibits complex flow structures. In addition to the pore-scale Reynolds number being sufficiently small, the condition for the existence of Darcy flow also requires that the characteristic scale of the flow field structure be much larger than the pore scale (Hewitt Reference Hewitt2020; De Paoli Reference De Paoli2023). It should be noted that the horizontal scale of the plumes is generally much smaller than their vertical scale. Therefore, in this section, we further examine the variations in the horizontal characteristic scale of the flow field.

Figure 5.

Time-averaged spectra of the mass concentration at $ d/2 \leqslant |z-0.5| \leqslant 2/D$ (excluding obstacles). The sampled data for fast Fourier transform (FFT) comes from the statistical steady state of the cases with $\phi =0.64$ . The triangles indicate the first dominant wavenumber, while the circles indicate the second dominant wavenumber. The wavenumber is re-scaled as $\hat {k}=k_x\Gamma /2\pi$ .

First, a simple method of finding the characteristic wavelength $\lambda _x$ , or the corresponding wavenumber $k_x$ , is the fast Fourier transform (FFT). Figure 5 shows the horizontal wavenumber spectra for all cases with $\phi =0.64$ . During sampling, we selected data from the middle height region of the domain to avoid the influence of the upper and lower boundaries. Since obstacles are arranged at $z=0.5$ , we selected the interval of $ d/2 \leqslant |z-0.5| \leqslant 2/D$ , which is a region within the height of one basic unit, excluding the obstacles. Additionally, for each case, we performed FFT on the instantaneous concentration field for at least 100 moments over the statistical time $t_s$ and then averaged them over time. The wavenumber is rescaled as $\hat {k}=k_x\Gamma /2\pi$ , thus it directly represents the number of periodic structures within the domain width. In figure 5, we use triangles and circles to denote the first and second dominant wavenumbers, respectively. From figure 5(a), we can see that for the smallest Rayleigh number ( ${\ {Ra}}=10^7$ ) considered here, the wavenumber $\hat {k}=4$ is absolutely dominant, corresponding to the four stable pairs of convection cells in figure 3(a). However, when $\ {Ra}$ increases to $4\times 10^7$ and $10^8$ , $\hat {k}=2$ replaces $\hat {k}=4$ as the dominant wavenumber, with their proportions being relatively close. In addition, other small wavenumbers, such as $\hat {k}=1$ and $\hat {k}=3$ , also have significant shares. This indicates that the flow field is tending towards a more chaotic state. When $\ {Ra}$ increases to $2\times 10^8$ , the first dominant wavenumber is $\hat {k}=1$ , corresponding to the highly chaotic and random flow state in figure 3(b), where there is no obvious periodic structure. For the three cases in figure 5(b), the shares of the small wavenumbers are relatively even, indicating that there is no stable periodic structure in the flow field at this stage, and plumes of multiple scales are distributed chaotically in both time and space. When $\ {Ra}$ increases further to $2\times 10^9$ , the wavenumber $\hat {k}=2$ becomes dominant and its proportion far exceeds that of other wavenumbers, as shown in figure 5(c). This corresponds to the two pairs of convection cells in figure 3(c). After this, as $\ {Ra}$ increases, the dominant wavenumber remains at $\hat {k}=2$ and the intensity of all modes gradually decreases. When $\ {Ra}$ reaches $10^{11}$ , the intensity of the dominant wavenumber is quite low, implying a weak disturbance in the concentration field, as shown in figure 3(d). Overall, as $\ {Ra}$ or ${\ {Ra}}_D$ increases, the dominant wavenumber of the flow field exhibits a relatively complex, non-monotonic change, with multiple wavenumbers interacting with each other. Previous numerical simulations under the Darcy model indicate $k_x\sim {\ {Ra}}_D^{0.4}$ (Hewitt et al. Reference Hewitt, Neufeld and Lister2014). In current pore-scale simulations, due to the interaction between large-scale convection cells and small-scale plumes, FFT alone cannot provide similar conclusions.

Figure 6.

(a) Horizontal auto-correlation function $R_x(\delta _x)$ of the mass concentration at $ d/2 \leqslant |z-0.5| \leqslant 2/D$ (excluding obstacles). The sampled data come from the statistical steady state of each case with $\phi =0.64$ . The dashed lines denote the quadratic fitting of $R_x$ within $R_x \geqslant 0.9$ for cases of ${\ {Ra}}=10^{7}$ and ${\ {Ra}}=10^{11}$ . (b) Horizontal characteristic wavelength $\lambda _x$ versus the Rayleigh number $\ {Ra}$ . $\lambda _x$ is four times the $\delta _x$ value at $R_x=0$ from the quadratic fitting curve shown in panel (a).

To uniformly examine the variation pattern of the characteristic scale, we calculated the horizontal auto-correlation coefficient of the concentration field within the middle height region, namely:

(3.1) \begin{equation} R_x(\delta _x)=\frac {\langle \overline {c(x,z,t)c(x+\delta _x,z,t)} \rangle _{z_{mid}}}{\langle \overline { c^2(x,z,t)}\rangle _{z_{mid}}}. \end{equation}

For $\phi =0.64$ and $0.15$ , $z_{mid}$ is $d/2 \leqslant |z-0.5| \leqslant D/2$ ; for $\phi =0.36$ , $z_{mid}$ is $|z-0.5| \leqslant l/2$ . Figure 6(a) plots $R_x(\delta _x)$ for all cases with $\phi =0.64$ . The colours denote different $\ {Ra}$ . It can be seen that when ${\ {Ra}}=10^7$ , $R_x$ rapidly decreases from 1 to below 0 as $\delta$ increases, without any inflection points. This corresponds to the periodic structure shown in figure 3(a). If we assume that it satisfies a sinusoidal model (Hewitt et al. Reference Hewitt, Neufeld and Lister2014), that is, $c=c_0(z,t)\sin(2\pi x/\lambda _x)$ , then we can easily find $R_x(\lambda _x/4)=0$ . This means that we can extract four times the $\delta _x$ value at $R_x=0$ as the horizontal wavelength. However, as $\ {Ra}$ increases, $R_x$ gradually deflects above 0. When $\ {Ra}$ reaches approximately $10^9$ , a clear inflection point appears at approximately $R_x=0.8$ , corresponding to the emergence of large-scale structures (see figure 3 c). Before the inflection point, $R_x$ still reflects the characteristics of small-scale structures to some extent. Therefore, we extracted the $R_x$ values between 0.9 and 1 and performed a quadratic fitting, as shown by the dashed lines in figure 6(a). The value of $\delta _x$ at $R_x=0$ of the fitted curve can be approximately considered as a quarter of the wavelength of the small-scale structure. Using this method, we calculated the horizontal wavelengths for all cases, with the results plotted in figure 6(b). For three sets of cases with different porosities, the variation of $\lambda _x$ with $\ {Ra}$ is similar. That is, as $\ {Ra}$ increases, $\lambda _x$ first decreases rapidly, then remains unchanged, and even slightly increases, before continuing to decrease. The intermediate plateau phase reflects the complex transition of the flow regime from Darcy flow to RB convection. The increase in scale with increasing $\ {Ra}$ might be due to dispersion effects, which is also found in experiments (Liang et al. Reference Liang, Wen, Hesse and DiCarlo2018).

Figure 7.

(a) Horizontal characteristic wavelength $\lambda _x$ re-scaled by $D$ versus the Rayleigh–Darcy number ${\ {Ra}}_D$ . The vertical dashed line denotes ${\ {Ra}}_D=4000$ . (b) $\lambda _x$ versus the ratio of the pore space $l$ and concentration BL thickness. The vertical dashed line denotes $l/\delta _c=6$ . Both panels are in log-log coordinates. The open coloured symbols denote the Darcy regime, the grey symbols denote the transition regime and the solid symbols denote the RB regime.

Previous pore-scale simulations (Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020) indicate that the scale of interior plumes in the Darcy regime increases with increasing $D$ , namely the scale of the basic unit. In figure 7(a), we show the variation of $\lambda _x/D$ with ${\ {Ra}}_D$ . It can be seen that when ${\ {Ra}}_D$ is small, the cases with different $\phi$ largely overlap; and as $\phi$ decreases, the wavenumber $k_x=2\pi /\lambda _x$ and ${\ {Ra}}_D$ gradually approach the $0.4$ scaling (Hewitt et al. Reference Hewitt, Neufeld and Lister2014; Gasow et al. Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020), as indicated by the black dashed line in the figure. When ${\ {Ra}}_D$ reaches approximately 4000, $\lambda _x/D$ starts to stabilise at approximately 2 and does not continue to decrease, implying that the width of a single plume ( $\lambda _x/2$ ) remains comparable to the scale of one basic unit. This indicates that when the characteristic scale of the flow field structure decreases to the pore scale, the Darcy flow assumption starts to fail. Over a considerable range of ${\ {Ra}}_D$ values thereafter (up to approximately ${\ {Ra}}_D=10^5$ ), although the intrinsic scale of the plumes decreases, the dispersion effect caused by the porous structure leads to its actual scale always increasing to the size of the basic unit. We define the flow state at this stage as being in the transition regime, indicated by grey symbols in figure 7. Before this stage, the variation in characteristic scale follows the features of the Darcy flow, so we refer to the flow state at this stage as being in the Darcy regime.

After the transition regime, $\lambda _x$ begins to decrease continuously with increasing ${\ {Ra}}_D$ . At this stage, the cases with different $\phi$ separate, indicating that the characteristic scale is no longer controlled by $D$ . Similarly, ${\ {Ra}}_D$ , as the control parameter that describes the Darcy regime, is no longer applicable at this stage. From figure 3, we can observe that the flow field enters a state similar to RB convection. Our previous work found that when the thickness of the thermal boundary layer (BL) is less than the pore space $l$ , the heat transport begins to follow the classic $0.3$ scaling in RB convection (Liu et al. Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020). Therefore, we also examine the relationship between the concentration BL thickness $\delta _c$ and $l$ in current cases. Here, dimensionless thickness is defined as $\delta _c=\delta _c^*/H^*=1/(2{\ {Sh}})$ . Figure 7(b) shows the relationship between $\lambda _x$ and $l/\delta _c$ . We find that the transition regime does not end at $l/\delta _c=1$ , but at approximately $l/\delta _c=6$ , as indicated by the grey dashed line in the figure. After this, the cases with different $\phi$ overlap again and satisfy the linear relationship $\lambda _x \sim \delta _c/l$ . We refer to the flow state that satisfies this relationship as being in the RB regime, indicated by solid symbols. In the next section, we will show that $\ {Sh}$ and $\ {Ra}$ indeed satisfy the 0.3 scaling of RB convection in this regime.

Figure 8.

Zoom-in plots of the instantaneous concentration fields near the bottom plate. The two cases have the same ${\ {Ra}}=10^{11}$ and $\delta _c=0.002$ , but with different porosities: (a) $\phi =0.64$ ; (b) $\phi =0.36$ . The plumes within the dashed ellipse undergo mechanical dispersion due to the obstacles.

It is natural that $\lambda _x$ and $\delta _c$ satisfy a linear relationship because the scale of the plume shedding from the boundary layer is similar to the BL thickness (van der Poel et al. 2015). However, why is $\lambda _x$ inversely related to $l$ ? Figure 8 shows the instantaneous concentration fields near the bottom boundary for two cases with $\phi =0.64$ and $0.36$ . These two cases have the same Rayleigh number ${\ {Ra}}=10^{11}$ and their Sherwood numbers $\ {Sh}$ are similar, resulting in almost the same BL thickness $\delta _c$ , both of which are smaller than the corresponding pore space $l$ . From figure 8(a), we can see that when $l$ is much larger than $\delta _c$ , the plumes can grow relatively freely. However, when they encounter obstacles, mechanical dispersion occurs, leading to an increase in the horizontal scale, as seen in the dashed ellipses in the figure. A naive idea is to assume that one plume maintains its volume when it encounters obstacles. If it has an intrinsic vertical scale $\lambda _z$ , we have $\lambda _z \delta _c=\lambda _x l$ , from which we can derive $\lambda _x \sim \delta _c/l$ . When $l$ and $\delta _c$ are comparative, as shown in figure 8(b), the constraint on the plumes by the obstacles becomes very strong. At this stage, the plumes can disperse to the scale of a basic unit and the flow field is still in the transition regime. This explains why the transition to the RB regime requires $l$ to be much larger than $\delta _c$ . Additionally, we should note that in the RB regime, although the porous structure no longer affects the BLs and, therefore, does not affect the transport efficiency, it still impacts the small-scale structure of the internal flow field.

4.1. Mass transport efficiency

Now, we examine the efficiency of mass transport in different regimes. In the Darcy model, classical theory posits that when flow enters the high Rayleigh regime (hereafter referred to as the high- ${\ {Ra}}_D$ Darcy regime), $\ {Sh}$ and ${\ {Ra}}_D$ satisfy a linear relationship, i.e. ${\ {Sh}} \sim {\ {Ra}}_D$ (Howard Reference Howard1966; Doering & Constantin Reference Doering and Constantin1998). Recent two-dimensional Darcy-scale simulations have corroborated this conclusion well (Hewitt et al. Reference Hewitt, Neufeld and Lister2014). However, both pore-scale simulations and experiments have observed nonlinear scaling. For example, Gasow et al.’s (Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020) pore-scale simulations found ${\ {Sh}} \sim {\ {Ra}}_D^{0.8}$ for $\phi =0.56$ and Liang et al.’s (Reference Liang, Wen, Hesse and DiCarlo2018) one-sided experiments found ${\ {Sh}} \sim {\ {Ra}}_D^{0.75}$ . This indicates that the actual pore-scale structures can influence mass transfer.

Figure 9.

Time-averaged Sherwood number $\overline {{\ {Sh}}}$ versus the Rayleigh–Darcy number ${\ {Ra}}_D$ . The cases with $\phi =0.15{-}0.64$ are from the current study, where the obstacles are impermeable to the species, while the cases with $\phi =0.75{-}0.92$ are from Liu et al. (Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020) and Xu et al. (Reference Xu, Xu and Xi2023), where the obstacles are permeable and impermeable to heat, respectively. For the current cases, the open coloured symbols denote the Darcy regime, the grey symbols denote the transition regime and the solid symbols denote the RB regime. The yellow area denotes the high- $Ra_D$ Darcy regime with $1300 \leqslant Ra_D \leqslant 4000$ ; the exponent $\gamma$ is calculated by fitting the data within this area, which equals 0.81, 0.92 and 0.97 for $\phi =0.64$ , 0.36 and 0.15, respectively.

In the current study, we similarly observed nonlinear relationship for $\ {Sh}$ and ${\ {Ra}}_D$ , as shown in figure 9. It can be seen that within the Darcy regime, the data approach the linear scaling (indicated by the black dashed line). However, we performed a linear fit on the data where $1300 \leqslant {\ {Ra}}_D \leqslant 4000$ , namely the high- ${\ {Ra}}_D$ Darcy regime, and found ${\ {Sh}} \sim {\ {Ra}}_D^{0.81}$ for $\phi =0.64$ , ${\ {Sh}} \sim {\ {Ra}}_D^{0.92}$ for $\phi =0.36$ and ${\ {Sh}} \sim {\ {Ra}}_D^{0.97}$ for $\phi =0.15$ . These are consistent with the conclusions of Gasow et al. (Reference Gasow, Lin, Zhang, Kuznetsov, Avila and Jin2020, Reference Gasow, Kuznetsov, Avila and Jin2021, Reference Gasow, Kuznetsov and Jin2022). However, with further increases in ${\ {Ra}}_D$ , our current cases reveal a transition from ${\ {Sh}} \sim {\ {Ra}}_D$ to ${\ {Sh}} \sim {\ {Ra}}^{0.3}$ (note that ${\ {Ra}}_D \sim {\ {Ra}}$ when $Da$ remains constant). This transition begins at approximately ${\ {Ra}}_D=4000$ , which is more easily observed in figure 7(a). Furthermore, within the transition regime, the scaling exponent of $\ {Sh}$ with ${\ {Ra}}_D$ gradually changes, making the initial phase of the transition regime potentially imperceptible. In figure 8 of Gasow et al. (Reference Gasow, Kuznetsov and Jin2022), it can also be observed that as ${\ {Ra}}_D$ approaches $10^5$ , there is a trend of decreasing $\ {Sh}$ . The underlying physical interpretation of the transition may be the influence of the porous structure on the flow structure. The transition starts when the characteristic scale of the flow structure is comparable to the size of the basic unit. When the scale of the plumes, or the thickness of the boundary layer, is much smaller than the size of the pore space $l$ , the effect of the porous media on mass transport can be neglected and $\ {Sh}$ with $\ {Ra}$ satisfies the 0.3 scaling of RB convection. For smaller $\phi$ (corresponding to smaller $l$ ), the Rayleigh number required to reach the RB regime is larger. The current cases with $\phi = 0.15$ have not yet reached the RB regime, which would require $\ {Ra}$ to exceed the order of $10^{13}$ , posing a significant challenge under current computational conditions. Finally, since thermal turbulence has long been a subject of interest among researchers and has been extensively explored in various structures (e.g. Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020; Wang et al. Reference Wang, Liu, Zhou and Sun2023, Reference Wang, Liu, Lai and Sun2024), we compare the results obtained from the previous work (Liu et al. Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020; Xu et al. Reference Xu, Xu and Xi2023) on thermal convection in porous media, as shown by the purple lines and orange symbols in the figure. The differences between the two primarily lie in two aspects: whether the obstacles are adiabatic and whether the obstacles are arranged in a regular pattern. In these cases, the Prandtl number ( $Pr=0.7{-}5.3$ ) is significantly smaller than the currently set Schmidt number ( ${\ {Sc}}=400$ ) and the porosity is also higher ( $0.75 \leqslant \phi \leqslant 0.92$ ), resulting in a much earlier transition. Consequently, the flow enters the RB regime at a relatively low Rayleigh–Darcy number. In the RB regime, the data for both types of thermal convection align closely, indicating that the adiabatic nature and regular arrangement of the obstacles have minimal impact on the overall thermal transport efficiency. However, when ${\ {Ra}}_D$ is very small, Xu et al.’s (Reference Xu, Xu and Xi2023) results are larger than those of Liu et al. (Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020). This is likely because irregularly arranged obstacles disrupt the original quasi-steady structure, leading to more intense flow. This phenomenon is also observed in our test cases, as shown in Appendix A.

Figure 10.

(a) Time-averaged Sherwood number $\overline {{\ {Sh}}}$ versus the Darcy number $Da$ . (b) $\overline {{\ {Sh}}}$ versus the Péclet number $Pe$ . (c) Ratio of the root-mean-square velocity $\boldsymbol {u}^*_{rms}$ to the characteristic velocity $U_K^*$ defined by the permeability versus $Pe$ . (d) Ratio of the effective diffusivity $\kappa _e$ to the molecular diffusivity $\kappa$ versus $Pe$ . The open coloured symbols denote the Darcy regime, the grey symbols denote the transition regime and the solid symbols denote the RB regime.

We further examine the effect of the Darcy number $Da$ itself on transport efficiency. In figure 10(a), we plot the variation of $\overline {{\ {Sh}}}$ with $Da$ under different $Ra$ . The results align with expectations: $\overline {{\ {Sh}}}$ decreases as $Da$ decreases, indicating that transport efficiency decreases with decreasing permeability. Moreover, a smaller $Da$ allows the flow to remain in the Darcy regime for a longer period. For example, at ${\ {Ra}}=10^{10}$ , the case with ${\ {Da}}=9.1 \times 10^{-8}$ remains in the Darcy regime, while the case with ${\ {Da}}=1.7 \times 10^{-6}$ has already transitioned to the transition regime and the case with ${\ {Da}}=2.6 \times 10^{-5}$ has entered the RB regime. As previously mentioned, the Darcy model is valid when the characteristic scale of the flow field is much larger than the pore scale. This condition is equivalent to the pore-scale Péclet number $Pe$ being much less than 1 (Hewitt Reference Hewitt2020; De Paoli Reference De Paoli2023). Here, $Pe$ is defined as

(4.1) \begin{equation} Pe=\frac {U^*_K\sqrt {K^*}}{\kappa }={\ {Ra}}_DDa^{1/2}, \end{equation}

in which $U^*_K=g\beta \varDelta _c^*K^*/\nu$ is the characteristic buoyancy velocity in the Darcy model. Note that $Pe$ represents the ratio of the pore scale $\sqrt {K^*}$ to the characteristic length scale $\kappa /U_K^*$ . Figure 10(b) illustrates the variation of $\overline {{\ {Sh}}}$ with $Pe$ . It can be observed that for cases in the Darcy regime, $Pe$ is not significantly less than 1 and can reach magnitudes of the order of $O(10)$ . However, the results in figure 7(a) indicate that for these cases, the horizontal length scale $\lambda _x^* \gtrsim 2D^* \gg \sqrt {K^*}$ . This apparent contradiction arises from the underestimation of the actual flow scale by the length scale $\kappa /U_K^*$ .

To clarify this point, we examine the discrepancies between the velocity $U_K^*$ and the diffusivity $\kappa$ compared with their actual effective values. First, we calculate the ratio of the root-mean-square velocity $\boldsymbol {u}^*_{rms}$ to $U_K^*$ , as shown in figure 10(c). It is evident that within the Darcy regime, this ratio remains relatively constant. For higher porosities $\phi =0.64$ and 0.36, $U_K^*$ slightly overestimates $\boldsymbol {u}^*_{rms}$ , while for the smallest porosity $\phi =0.15$ , $U_K^*$ is nearly equal to $\boldsymbol {u}^*_{rms}$ . As the flow enters the transition regime, $U_K^*$ increasingly overestimates the actual velocity, signalling the breakdown of the Darcy model. Second, we define an effective diffusivity based on the actual velocity $\boldsymbol {u}^*_{rms}$ and characteristic length $\lambda _x^*$ , i.e.

(4.2) \begin{equation} \kappa _e=\lambda _x^*\boldsymbol {u}^*_{rms}. \end{equation}

We compute the ratio of $\kappa _e$ to $\kappa$ , as shown in figure 10(d). The results reveal that, for all cases within the Darcy regime, the magnitude of $\kappa _e/\kappa$ lies in the range of $O(10)$ – $O(100)$ . This indicates that in pore-scale simulations, the effective diffusivity is generally much larger than the molecular diffusivity, which can be essentially attributed to the effects of mechanical dispersion. This phenomenon is the primary reason why the length scale in simulations exceeds expectations. Note that in the transitional and RB regimes, the dispersion effect is further amplified, implying that the length scale remains influenced by the porous structure. This observation is consistent with the conclusions in § 3.2.

4.2. Energy dissipation rate

Figure 11.

Normalised dissipation rates for (a) concentration and (b) kinetic energy versus the distance from the wall. The distance is re-scaled by the size $D$ of the basic unit. The dissipation rates are averaged both in time and space ( $(n-1)D\leqslant z \leqslant nD, n=1,2,3\ldots$ ). The open coloured symbols denote the cases with $Ra_D=2.6 \times 10^3$ ( $\phi =0.64$ and 0.36) and $Ra_D=2.7 \times 10^3$ ( $\phi =0.15$ ); the grey symbols denote the cases with $Ra_D=2.6 \times 10^4$ ( $\phi =0.64$ ), $Ra_D=3.4 \times 10^4$ ( $\phi =0.36$ ) and $Ra_D=3.6 \times 10^4$ ( $\phi =0.15$ ); the solid symbols denote the cases with $Ra_D=2.6 \times 10^6$ ( $\phi =0.64$ ) and $Ra_D=1.7 \times 10^6$ ( $\phi =0.36$ ).

We further investigate the energy dissipation in different regimes. The dimensional energy dissipation rates for concentration and kinetic energy are defined as follows:

(4.3) \begin{equation} \epsilon ^*_c=\kappa \sum _i \left [ \frac {\partial c^*}{\partial x_i^*}\right ]^2, \quad \epsilon _u^*=\frac {1}{2}\nu \sum _{ij} \left [ \frac {\partial u_j^*}{\partial x_i^*}+\frac {\partial u_i^*}{\partial x_j^*}\right ]^2. \end{equation}

The corresponding dimensionless forms read

(4.4) \begin{equation} \epsilon _c=\frac {1}{\sqrt {RaSc}}\sum _i \left [ \frac {\partial c}{\partial x_i}\right ]^2, \quad \epsilon _u=\frac {1}{2}\sqrt {\frac {Sc}{Ra}}\sum _{ij} \left [ \frac {\partial u_j}{\partial x_i}+\frac {\partial u_i}{\partial x_j}\right ]^2. \end{equation}

In figure 11, we plot the vertical distribution of the two energy dissipation rates. Here, the energy dissipation rates are time-averaged and spatially averaged over vertical intervals of $D$ , namely the size of the basic unit. Additionally, the vertical coordinate is also rescaled using $D$ . For each regime and porosity, we selected cases with similar ${\ {Ra}}_D$ values. In figure 11(a), it can be observed that, after normalising by the value at the first point of each group (denoted by superscript 1), the distribution of $\epsilon _c$ for cases with similar $Ra_D$ is consistent. The dissipation rate rapidly decreases with increasing distance from the wall, and this decrease becomes more pronounced as $Ra_D$ increases. In the Darcy regime, $\epsilon _c$ decreases to approximately $4\,\%$ of $\epsilon _c^1$ at a distance of three basic units from the wall. For the transition regime and RB regime, this ratio is approximately $2\,\%$ and $0.5\,\%$ , respectively. Therefore, concentration energy dissipation is concentrated mainly within the distance $3D$ from the wall. Gasow et al. (Reference Gasow, Kuznetsov and Jin2022) also reported that in their numerical simulations, energy dissipation within 5 basic units of the wall accounted for $93\,\%$ of the total dissipation. The distribution of kinetic energy dissipation, however, shows a different pattern. As shown in figure 11(b), at lower ${\ {Ra}}_D$ values, there is a slight increase in $\epsilon _u$ near the centre, reaching approximately $18\,\%$ of $\epsilon _u^1$ . This is caused by the chaotic plumes observed in the high- $Ra_D$ Darcy regime and the transition regime. Numerous plume structures in the bulk area cause significant kinetic energy dissipation at the boundaries of obstacles. Here, $\epsilon _c$ is less affected because the normal concentration gradient at the obstacles’ boundaries is zero. Moreover, for different porosities, the change in $\epsilon _u$ with $Ra_D$ exhibits a lag effect. For larger porosities, as $Ra_D$ increases, $\epsilon _u$ increases later in the central region before eventually decreasing. In the RB regime, $\epsilon _u$ continues to decrease with increasing distance from the wall. At the fifth basic unit position, it has dropped to approximately $1\,\%$ of $\epsilon _u^1$ .

Figure 12.

Typical snapshots of (a,c,e) concentration energy dissipation rate log $_{10}\epsilon _c$ and (b,d,f) kinetic energy dissipation rate log $_{10}\epsilon _u$ for the cases with $\phi =0.64$ . The yellow dashed lines in panels (a,c,e) denote the concentration BL.

The distribution of dissipation rates throughout the entire flow field allows us to see the differences between the various regimes more clearly, as shown in figure 12. Here, we selected three cases represented by circular symbols ( $\phi =0.64$ ) from figure 11. In the high- $Ra_D$ Darcy regime (figures 12 a and 12 b), the chaotic plumes lead to an increase in both dissipation rates. For the concentration, dissipation is more pronounced near the boundaries and extends several basic units into the bulk region. In contrast, significant kinetic energy dissipation shows a larger region throughout the field, with most of the contribution coming from plume structures interacting with the boundaries of obstacles. In the transition regime (figures 12 c and 12 d), plume scales become smaller and large-scale structures gradually form. This causes concentration energy dissipation to be more concentrated near the boundaries, with very low dissipation across a wide area in the bulk region. In contrast, for kinetic energy dissipation, the more intense flow in the bulk region leads to increased dissipation, even exceeding that near the boundaries. As the flow enters the RB regime (figures 12 e and 12 f), stable large-scale convective rolls result in an even larger low- $\epsilon _c$ area within the bulk region, with dissipation concentrated mainly along the edges of the convective rolls. At this stage, kinetic energy dissipation throughout the field further increases and contributions near the boundaries become more significant because of the presence of chaotic small-scale plumes near the boundaries.

In the next section, we will analyse the $Sh {-} Ra$ relationship based on the Grossmann and Lohse (GL) theory (Grossmann & Lohse Reference Grossmann and Lohse2000, 2001), which primarily focuses on dividing energy dissipation rates into BL contributions and bulk contributions. Therefore, here, we first examine the BL thickness. For the concentration field, we calculated the dimensionless BL thickness using $\delta _c=1/(2{\ {Sh}})$ . In figure 12(a,c,e), concentration boundary layers (BLs) are marked with yellow dashed lines. It can be observed that, in the Darcy regime, $\delta _c$ is approximately 0.028, exceeding the pore width $l$ , but smaller than the basic unit size $D$ . In this case, the dissipation rate within BLs is relatively large. In the transition regime, $\delta _c$ is reduced to less than $l$ , approximately 0.008. In the RB regime, $\delta _c$ becomes very small, approximately 0.004, resulting in a minimal contribution to the overall dissipation. For the kinetic BL, due to the no-slip condition at the obstacles, the BL thickness at both the upper and lower plates is always less than $l$ . In our subsequent theoretical analysis, we consider that BLs exist at the boundaries of all obstacles throughout the field. The definitions of the bulk region and the BL region differ from the traditional ones in this context. Specific details will be introduced in the next section.