From Darcy convection to free-fluid convection: pore-scale study of density-driven currents in porous media

Abstract

We conducted a series of pore-scale numerical simulations on convective flow in porous media, with a fixed Schmidt number of 400 and a wide range of Rayleigh numbers. The porous media are modeled using regularly arranged square obstacles in a Rayleigh–Bénard (RB) system. As the Rayleigh number increases, the flow transitions from a Darcy-type regime to an RB-type regime, with the corresponding $Sh$–$Ra_D$ relationship shifting from sublinear scaling to the classical 0.3 scaling of RB convection. Here, $Sh$ and $Ra_D$ represent the Sherwood number and the Rayleigh–Darcy number, respectively. For different porosities, the transition begins at approximately $Ra_D = 4000$, at which point the characteristic horizontal scale of the flow field is comparable to the size of a single obstacle unit. When the thickness of the concentration boundary layer is less than approximately one-sixth of the pore spacing, the flow finally enters the RB regime. In the Darcy regime, the scaling exponent of $Sh$ and $Ra_D$ decreases as porosity increases. Based on the Grossman–Lohse theory (J. Fluid Mech. vol. 407, 2000, pp. 27–56; Phys. Rev. Lett. vol. 86, 2001, p. 3316), we provide an explanation for the scaling laws in each regime and highlight the significant impact of mechanical dispersion effects during the development of the plumes. Our findings provide some new insights into the validity range of the Darcy model.

1. Introduction

Porous media refer to materials composed of numerous frameworks that create many tiny voids. When a fluid fills these voids, convection may occur under the influence of gravity due to density differences. This density-driven convection in porous media exhibits properties that differ from those of free fluid convection and is widely present in nature and engineering applications, such as the formation of sea ice and the oil recovery from geological formations (Farajzadeh et al. 2012; Miah et al. 2018; Anderson et al. 2022). In recent years, due to the rise of carbon dioxide (CO $_2$ ) geological sequestration, convection in porous media has become a research hotspot (Huppert & Neufeld 2014; Bachu 2015; De Paoli 2021). This technology involves capturing CO $_2$ generated from industrial processes and burying it underground to reduce carbon emissions. When CO $_2$ is injected into deep saline aquifers, it enters a supercritical state similar to a fluid and initially rises to an impermeable layer due to its lower density. Subsequently, CO $_2$ diffuses horizontally and gradually dissolves in the underlying saline water. This dissolved liquid is denser than the surrounding fluid, thus driving convection. The dissolution and convection process of CO $_2$ in saline aquifers is crucial for preventing leakage and ensuring its long-term storage. Therefore, understanding the mass transport efficiency in porous media convection has become a key focus of research (Hewitt 2020; De Paoli 2023).

The commonly used model for convective flow in porous media is based on Darcy’s law (hereafter referred to as the Darcy model) (Nield & Bejan 2017). The smallest analytical flow domain in the Darcy model includes multiple pores and is known as the representative elementary volume (REV). The REV scale is generally much larger than the pore scale, but much smaller than the characteristic flow scale. Consequently, the Darcy model solves for macroscopic variables that are volume-averaged, neglecting information at the pore scale. In this model, the most important governing parameter is the Rayleigh–Darcy number ${\ {Ra}}_D$ . Early theoretical studies suggest that for asymptotically large ${\ {Ra}}_D$ , the Sherwood number $\ {Sh}$ , which represents the efficiency of mass transport, scales linearly with ${\ {Ra}}_D$ (Malkus 1954; Howard 1966). With the improvement of computational capabilities in recent years, numerous numerical simulations for the high- $Ra_D$ Darcy model have emerged (Hewitt et al. 2012, 2013, 2014; Pirozzoli et al. 2021; De Paoli et al. 2022). The results of two-dimensional (2-D) simulations confirm the linear scaling of $\ {Sh}$ versus ${\ {Ra}}_D$ . However, three-dimensional (3-D) simulations have revealed an additional sublinear term. This discrepancy is primarily due to the fact that in 3-D simulations, ${\ {Ra}}_D$ has not reached a sufficiently high value for the flow to enter the ultimate regime (De Paoli 2023). Recently, Zhu et al. (2024) successfully used the Grossman–Lohse (GL) theory (Grossmann & Lohse 2000, 2001) to explain the differences between 2-D and 3-D results.

Although numerical simulations align with the theory, numerous laboratory experiments have consistently found that the scaling exponent of $Sh$ versus $Ra_D$ is always less than 1 (Neufeld et al. 2010; Backhaus et al. 2011; Wang et al. 2016; Liang et al. 2018). This indicates that pore-scale effects cannot be neglected when considering mass transport in convective flows through actual porous media. To incorporate these effects into simulations, there are generally two approaches. One is to introduce additional terms into the Darcy model to account for the interactions between solid structures and the fluids. For example, when the flow velocity within the pores is relatively high, the drag exerted by solid obstacles on the flow becomes non-negligible, requiring consideration of the Forchheimer term (Joseph et al. 1982; Nield & Bejan 2017; Jin & Kuznetsov 2024). In recent years, there has been increased attention to mechanical dispersion, which refers to the alteration of flow direction and further mixing of solutes due to pore structures (Saffman 1959; Dentz et al. 2023). This dispersion effect produces results similar to those caused by molecular diffusion; therefore, the two are often considered together under the term hydrodynamic dispersion (De Paoli 2023). In the Darcy simulations, to account for the effects of dispersion, the molecular diffusion coefficient is typically replaced by the Fickian dispersion tensor (Bear 1961; Hidalgo & Carrera 2009; Ghesmat et al. 2011). Recent simulation (Wen et al. 2018) reported a fan-flow structure caused by dispersion effects, which indeed influence mass transport. However, because these models introduce numerous parameters, extensive experiments and numerical simulations are still required to calibrate and validate them.

Another approach involves conducting pore-scale simulations, which solve the original Navier–Stokes (NS) equations within the pores, rather than using the Darcy model. This places extremely high demands on computational capabilities, particularly for cases with low porosity. In our previous work (Liu et al. 2020, 2021), we conducted 2-D pore-scale simulations on thermal convection with relatively high porosities and found that heat transfer first increases and then decreases as porosity decreases. This is due to the fact that obstacles simultaneously hinder convective flow and improve the continuity of the flow structure. Xu et al.’s (2023) simulations using impermeable obstacles also led to similar conclusions, although the plume dynamics were less coherent compared with those within permeable obstacles. Recently, Gasow et al.’s (2020, 2021, 2022) conducted a series of 2-D pore-scale simulations at large Schmidt numbers with porosities as low as 0.09, discovering that the scaling exponent of $Sh$ versus $Ra_D$ decreases with increasing porosity, which aligns with experimental results. Moreover, they pointed out that the discrepancy between pore-scale simulations and Darcy simulations might be due to the latter’s omission of the momentum dispersion term.

One concern about pore-scale simulations is that as the control parameter, i.e. the Rayleigh number $\ {Ra}$ , increases, the characteristic scale of the flow field gradually decreases and the flow may eventually transition into a state similar to Rayleigh–Bénard (RB) convection that is no longer influenced by the porous media. Both simulations (Liu et al. 2020; Xu et al. 2023) and experiments (Ataei–Dadavi et al. 2019) have shown that under extremely high $\ {Ra}$ , the heat transport efficiency follows the classical 0.3 scaling of RB convection (Grossmann & Lohse 2000, 2001). For solutes with much lower diffusivity, would the flow eventually transition to an RB state as well? Gasow et al.’s (2020, 2021, 2022) simulations did not reveal any clear indications of this, possibly because the control parameter was not large enough. In this paper, we conduct systematic 2-D pore-scale simulations over a wide range of $\ {Ra}$ to investigate this issue, aiming to illustrate when the flow transitions from a Darcy-type to an RB-type regime and to thoroughly examine the characteristics of the transitional phase. We also aim to use the GL theory to analyse the transport scaling laws in different regimes. It should be noted that our pore-scale analysis fundamentally differs from the work of Zhu et al. (2024), as their study was conducted on the macroscopic (Darcy) scale.

This paper is organised as follows. In § 2, we introduce the governing equations and the numerical settings. Then, in §§ 3 and 4, we present the numerical results, including the variations of flow structures and mass transfer, respectively. In § 5, we provide a theoretical analysis based on the GL theory. Finally, we give the conclusions in § 6.

2. Problem formulation

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 2017). 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 2017):

(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 2020).

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 2017):

(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 2012). 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. (2020), 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 2024). 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. (2020, 2021, 2022), 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 2005). 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. 2017; Vanella & Balaras 2020), and was also employed in our previous simulations involving circular obstacles (Liu et al. 2020). 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. 2015, 2022; Li & Yang 2023). Basically, the finite difference method and the fractional time-step method are employed (Ostilla-Mónico et al. 2015). 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. 2022). In the current study, we fix ${\ {Sc}}=400$ , a value close to the practical conditions of CO $_2$ sequestration (Huppert & Neufeld 2014). 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 2017). 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. 2012; Wagner & Shishkina 2013; Shishkina 2021). 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 1983), 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. 2010). 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. Flow structure

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. 2014; Nield & Bejan 2017; Gasow et al. 2020): 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. 2018), 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 2024). 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 2020; De Paoli 2023). 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 2020; De Paoli 2023). 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. 2014). 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. 2014), 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. 2018).

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. 2020) 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. 2014; Gasow et al. 2020), 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. 2020). 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. Mass transfer and energy dissipation

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 1966; Doering & Constantin 1998). Recent two-dimensional Darcy-scale simulations have corroborated this conclusion well (Hewitt et al. 2014). However, both pore-scale simulations and experiments have observed nonlinear scaling. For example, Gasow et al.’s (2020) pore-scale simulations found ${\ {Sh}} \sim {\ {Ra}}_D^{0.8}$ for $\phi =0.56$ and Liang et al.’s (2018) 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. (2020) and Xu et al. (2023), 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. (2020, 2021, 2022). 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. (2022), 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. 2020; Wang et al. 2023, 2024), we compare the results obtained from the previous work (Liu et al. 2020; Xu et al. 2023) 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 (2023) results are larger than those of Liu et al. (2020). 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 2020; De Paoli 2023). 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. (2022) 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 2000, 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.

5. Application of GL theory

The key idea of GL theory is to decompose the dissipation rates into contributions from the BLs and the bulk region, i.e.

(5.1) \begin{equation} \epsilon _u^*=\epsilon _{u,BL}^*+\epsilon _{u,bulk}^*,\quad \epsilon _c^*=\epsilon _{c,BL}^*+\epsilon _{c,bulk}^*. \end{equation}

Depending on the regime, the dissipation rates may be dominated by the BL part or by the bulk part, leading to the corresponding estimation formulae. Moreover, in a closed RB system, there exist exact relations (Grossmann & Lohse 2000, 2001):

(5.2) \begin{equation} \langle \overline {\epsilon _u^*} \rangle _V=\frac {\nu ^3}{H^{*4}}({\ {Sh}}-1){\ {Ra}}{\ {Sc}}^{-2},\quad \langle \overline {\epsilon _c^*} \rangle _V=\kappa \frac {\varDelta _c^*}{H^{*2}}{\ {Sh}}. \end{equation}

Here, the overline represents time-averaging and the bracket $\langle \cdot \rangle _V$ indicates averaging over the entire domain. It should be noted that the presence of obstacles does not affect the validity of the exact relations (5.2).

In our previous work on thermal convection in porous media (Liu et al. 2020), we found that additional laminar-type BLs exist in the channels between obstacles, leading to the conclusion that kinetic energy dissipation is dominated by the numerous obstacle BLs throughout the domain. By combining the corresponding estimation formula with the exact relation (5.2), we obtain

(5.3) \begin{equation} {\ {Sh}} \approx \alpha \cdot \phi \left(\frac {H^*}{l^*}\right)^2{\ {Sc}}^2\ {Re}^2{\ {Ra}}^{-1}+1. \end{equation}

Here, $\alpha$ is an empirical coefficient. Note that we have replaced the Nusselt number $Nu$ and the Prandtl number $Pr$ in the original formula ((5.3) of Liu et al. 2020) with the Sherwood number $\ {Sh}$ and Schmidt number $\ {Sc}$ , respectively, and the Reynolds number $\ {Re}$ is defined based on the domain height $H^*$ (2.5). More details can be found in that paper. Although this formula was derived for thermal convection, it theoretically also holds for the current cases with large $\ {Sc}$ . In figure 13, we plot the relationship between $\ {Sh}$ and $\ {Ra}$ , along with (5.3) represented by the solid lines. Here, we set the coefficient $\alpha$ to 10 (compared with 8 by Liu et al. 2020). It can be observed that (5.3) matches the DNS results remarkably well across almost all regimes, including the high- ${\ {Ra}}_D$ Darcy regime, transition regime and RB regime. Together with previous work, (5.3) demonstrates high universality and is valid over a wide range of $\ {Sc}$ , $\ {Ra}$ and $\phi$ .

Figure 13. Time-averaged Sherwood number $\overline {{\ {Sh}}}$ versus the Rayleigh number $\ {Ra}$ . The open coloured symbols denote the Darcy regime, the grey symbols denote the transition regime and the solid symbols denote the RB regime. The solid lines denote (5.3) with $\alpha =10$ .

Nevertheless, (5.3) does not directly provide the relationship between $\ {Sh}$ and $\ {Ra}$ , as it also includes the Reynolds number $\ {Re}$ . Therefore, we need to further derive a new equation based on concentration dissipation. For large $\ {Sc}$ , as $\ {Ra}$ increases, the flow transitions from the so-called $I_\infty$ regime to $III_\infty$ regime (Grossmann & Lohse 2001), where $\epsilon _c$ shifts from being BL-dominated to bulk-dominated, consistent with our previous observations (figure 12). In the high- ${\ {Ra}}_D$ Darcy regime, $\epsilon _c$ is dominated by BLs and the governing equation for concentration within BLs can be simplified from (2.1b ) as follows:

(5.4) \begin{equation} u_x^*\partial _x c^*+u_z^*\partial _z c^*=\kappa \partial _z^2 c^*. \end{equation}

The magnitudes of the two terms on the left-hand side of the above equation are comparable due to the incompressibility condition (2.1c ) (Grossmann & Lohse 2000). However, in contrast to the classical RB configuration, the presence of the porous medium leads to the formation of numerous proto-plumes within the boundary layer, resulting in the magnitudes of $u_x^*$ and $u_z^*$ being comparable. This is discussed in detail in Appendix B. In summary, we can compute the order of magnitude of either term on the left-hand side of (5.4) to represent the overall magnitude. By using the balance between $u_z^*\partial _z c^*$ and $\kappa \partial _z^2 c^*$ and considering $\partial _z \sim 1/\delta _c^*$ , we can obtain

(5.5) \begin{equation} \frac {u_z^*}{\delta _c^*}\sim \frac {\kappa }{\delta _c^{*2}}. \end{equation}

If we define a vertical Reynolds number as $\ {Re}_z=\langle u^*_z \rangle ^{rms}_{BL} H^*/\nu$ , where $\langle \cdot \rangle ^{rms}_{BL}$ denotes the root-mean-square value calculated over the concentration BL region (including the obstacles), and consider ${\ {Sh}}=H^*/(2\delta _c^*)$ , we can conclude from (5.5) that

(5.6) \begin{equation} {\ {Sh}} \sim \ {Re}_z {\ {Sc}}. \end{equation}

In Appendix B, we demonstrate that $\ {Sh}$ and $\ {Re}_z$ indeed follow a linear scaling in the high- ${\ {Ra}}_D$ Darcy regime for various $\phi$ .

Figure 14. Time-averaged Sherwood number $\overline {{\ {Sh}}}$ versus (a) the time-averaged Reynolds number $\overline {Re}$ and (b) the Rayleigh number $\ {Ra}$ . The open coloured symbols denote the Darcy regime, the grey symbols denote the transition regime and the solid symbols denote the RB regime. The exponent $\zeta$ in panel (a) is calculated by fitting the DNS data, which equals 0.90, 0.96 and 0.98 for $\phi =0.64$ , 0.36 and 0.15, respectively. In panel (b), only cases in the high- $Ra_D$ Darcy regime are presented; the dashed line is calculated using (5.7) with $\zeta$ derived from panel (a); the solid line is given by Gasow et al.’s (2022) (15), where the geometric parameter $a$ takes values of 0.0125, 0.01 and 0.0092 for $\phi =0.64$ , 0.36 and 0.15, respectively.

If $\ {Re}_z \sim \ {Re}$ , then by combining (5.3) and (5.6), and considering ${\ {Sh}} \gg 1$ , we can obtain ${\ {Sh}} \sim {\ {Ra}}$ , consistent with the classical theory (Howard 1966). However, the vertival velocity within BLs originates from the proto-plumes. These plumes enter the bulk region and evolve into mega-plumes. During this process, the mechanical dispersion influences the changes in vertical velocity. Consequently, $\ {Re}_z \sim \ {Re}$ does not always hold when the porous structure has an impact on it. In figure 14(a), we plot the relationship between $Sh$ and $Re$ . By fitting the data points from the high- $Ra_D$ Darcy regime, the corresponding $Sh \sim Re^\zeta$ scaling is obtained, with exponents of 0.90, 0.96 and 0.98 for the three porosities. This suggests a nonlinear scaling $Re_z \sim Re^\zeta$ . In the RB regime, the data approach the $Sh \sim Re^{0.5}$ scaling, which can be derived from the classical GL theory in the $III_\infty$ regime (Grossmann & Lohse 2001). Substituting these two scaling laws into (5.3) and considering ${\ {Sh}} \gg 1$ , we can obtain

(5.7) \begin{equation} {\ {Sh}} \sim {\ {Ra}}^\gamma , \quad \gamma =\left \{ \begin{aligned} &\frac {\zeta }{2-\zeta } & \quad \textrm{high-}Ra_D \ \textrm {Darcy regime},\\ & \quad \frac {1}{3} & \quad \textrm{ RB regime}. \end{aligned} \right . \end{equation}

Here, $\zeta$ is an empirical coefficient related to the mechanical dispersion effect, reflecting the correlation between the vertical velocity within the BLs and the global velocity. As porosity decreases, $\zeta$ gradually approaches 1. Substituting its specific values into (5.7), we obtain the $\gamma$ values for the three porosities as 0.82, 0.92 and 0.96, respectively, which is consistent with the previous results (see figure 9). It should be noted that Gasow et al. (2022) also proposed an empirical correlation, i.e. ${\ {Sh}}=a{\ {Ra}}_D^{1-0.2\phi ^2}+1$ , where $a$ is a pore-scale geometric parameter. In figure 14(b), we plot ${\ {Sh}}{-}{\ {Ra}}$ for the cases in the high- ${\ {Ra}}_D$ regime, meanwhile, we show the results of (5.7) (denoted by dashed lines) and Gasow et al. (2022) (denoted by solid lines). It can be seen that both formulae agree well with the simulation results, while the exponents $\gamma$ predicted by Gasow et al. (2022) are slightly higher. Moreover, the range of the fitted coefficient $a$ is consistent with the interval obtained by Gasow et al. (2022). Overall, both formulae suggest that the pore-scale structure influences mass transport. What is particularly notable about the current study is that we successfully extend the GL theory to porous medium convection and provide a physical explanation for the nonlinear scaling of mass transport.

6. Conclusion

In summary, we conducted a series of 2-D pore-scale simulations to investigate density-driven convection in porous media. The porous structure was modelled by arranging obstacles in a regular pattern within the RB region, with porosities $\phi$ of 0.64, 0.36 and 0.15. With the Schmidt number set to 400, owing to the enhanced capabilities of the dual-resolution technique, we were able to set the values of the Rayleigh–Darcy number ${\ {Ra}}_D$ over a broad range, form $10^2$ up to $10^6$ . As ${\ {Ra}}_D$ increases, the flow regime sequentially transitions through the Darcy regime, transition regime and RB regime. Correspondingly, the $Sh\sim Ra_D^\gamma$ relationship transitions from a sublinear scaling to the classic 0.3 scaling seen in RB convection. The exponent $\gamma$ in the high- ${\ {Ra}}_D$ Darcy regime increases with decreasing $\phi$ , i.e. $\gamma =0.81$ , 0.92 and 0.97 for $\phi =0.64$ , 0.36 and 0.15, respectively. A detailed analysis of the flow field structure revealed that when the plume width in the bulk region decreases to a scale comparable to the basic unit size $D$ , the pore-scale structure begins to influence the flow and the Darcy model starts to fail. The corresponding critical $Ra_D$ is approximately 4000 for all porosities. The flow finally enters the RB regime when the concentration BL thickness $\delta _c$ is less than $1/6$ of the pore space $l$ , meaning that the plume width is much smaller than the spacing between obstacles. At this stage, the porous structure has minimal impact on mass transport. Using GL theory, we successfully explained the changes in the scaling of ${\ {Sh}} \sim {\ {Ra}}_D^\gamma$ across different regimes and obtained a unified formula (5.7), which shows good agreement with the DNS results. We found that in the high- ${\ {Ra}}_D$ Darcy regime, the variation in $\gamma$ may originate from the mechanical dispersion, which causes velocity changes during plume development. Our findings shed some new light on the nonlinear scaling that emerges in mass transport when considering pore-scale effects.

In the Darcy regime, as ${\ {Ra}}_D$ increases, the concentration field transitions from large-scale convection rolls to chaotic small-scale plumes. At this stage, the pore-scale Reynolds number remains much smaller than 1, indicating a ‘pseudo-turbulent’ state. As ${\ {Ra}}_D$ further increases, the flow enters the transition regime, where large-scale periodic structures re-emerge, accompanied by fan-shaped structures caused by the dispersion effect. When ${\ {Ra}}_D$ becomes sufficiently high, the flow enters the RB regime and the entire flow field consists of both very small-scale plumes and large-scale convection rolls. By calculating the horizontal auto-correlation coefficient of the concentration field, we can systematically examine the variation in the horizontal characteristic length scale of small-scale structures. In the Darcy regime, when the horizontal wavelength $\lambda _x$ is re-scaled by $D$ , the cases with different $\phi$ collapse and the wavenumber approximately follows $k_x \sim {\ {Ra}}_D^{0.4}$ . In the transition regime, $\lambda _x$ remains roughly constant at $2D$ , indicating that individual plumes tend to remain within one basic unit. In the RB regime, although the porous structure has little effect on mass transport, $\lambda _x$ becomes proportional to $\delta _c$ and inversely proportional to $l$ . This is primarily due to the horizontal dispersion of the plumes after encountering obstacles. Further investigation revealed that mechanical dispersion affects all regimes, resulting in the actual characteristic scale $\lambda _x$ being larger than the estimated value $\kappa /U_K^*$ . The above results indicate that, in porous media convection, the basic unit size $D$ is a crucial intrinsic length scale. Similarly, when the vertical coordinate $z$ is re-scaled by $D$ , the vertical profiles of the concentration dissipation rate $\epsilon _c$ exhibit consistency at the same $Ra_D$ . However, the kinetic energy dissipation rate $\epsilon _u$ shows a certain degree of lag, mainly due to the additional kinetic energy dissipation caused by the boundaries of internal obstacles.

The starting point of our current work is geological CO $_2$ sequestration, where understanding the settling rate of CO $_2$ in saline aquifers is of significant practical importance. However, the current simulations still have some discrepancies with real-world conditions, including the non-random distribution of obstacles, the flow being confined to 2-D, and the fixed salinity distribution on both the upper and lower boundaries (where the real scenario is closer to one-sided convection). Despite these limitations, our current configuration provides a statistically steady state, allowing us to draw some meaningful conclusions. This is helpful for understanding convection in porous media itself, which can ultimately be applied to CO $_2$ sequestration. Our study also provides valuable insights into thermal convection in porous media. Through comparisons presented in this paper, it can be observed that, regardless of whether the obstacles are adiabatic, heat transfer tends to transition into the RB regime earlier and the breakdown of the Darcy model occurs at relatively low Rayleigh–Darcy number. This is likely due to the higher thermal diffusivity. We are also working on the case where obstacles are randomly distributed in a 3-D setting.