3.1. The global heat transport at $ \textit{Da}=1.21\times 10^{-5}$

Figure 2.

Results of the global heat transport. (a) Fluid Nusselt number $ \textit{Nu}_{\kern-1pt f}$ versus fluid Rayleigh number $ \textit{Ra}_{\kern-1pt f}$ . (b) Compensated plots of $ \textit{Nu}_{\kern-1pt f}/Ra_f^{0.3}$ versus $ \textit{Ra}_{\kern-1pt f}$ . (c) Medium Nusselt number $ \textit{Nu}_m$ versus Rayleigh–Darcy number $ \textit{Ra}$ . (d) Compensated plot of $ \textit{Nu}_m/Ra$ versus $ \textit{Ra}$ .

We first present the experimentally measured global heat transport data at a fixed Darcy number of $ \textit{Da} = 1.21 \times 10^{-5}$ using red diamonds in figure 2. The data span over four decades of Rayleigh–Darcy number, enabling us to capture the continuous and complex evolution of heat transport behaviour in PMC – from the conduction regime to the RBC-like regime. For ease of reference, literature data (including laboratory experiments, DOB simulations and pore-scale DNS) are also plotted with open symbols.

Figure 2(a) shows the relationships between the fluid Nusselt number $ \textit{Nu}_{\kern-1pt f}$ and the fluid Rayleigh numbers $ \textit{Ra}_{\kern-1pt f}$ , together with data measured in classical RBC (black squares) for comparison. Figure 2(b) shows the corresponding compensated plots of $ \textit{Nu}_{\kern-1pt f}/Ra_f^{0.3}$ . Firstly, we note that the global heat transfer data obtained from cells with different aspect ratios collapse excellently onto a single master curve (for both PMC and RBC, see table 1). This indicates that within the parameter range explored by the present study ( $\varGamma _{xz}\ge 0.29$ ), the geometrical confinement effect has a negligible influence on the global heat transport, which is consistent with classical RBC (Huang et al. Reference Huang, Kaczorowski, Ni and Xia2013; Xia et al. Reference Xia, Huang, Xie and Zhang2023). In particular, our data, previous experimental results and pore-scale DNS datasets all show a tendency to converge to a regime that features a single power law for sufficiently large $ \textit{Ra}_{\kern-1pt f}$ . Specifically, the scaling exponent approaches $\beta \sim 0.3$ , which is consistent with that of classical RBC (see the solid black symbols).

Figures 2(c) and 2(d) plot the same datasets in terms of the medium Nusselt number $ \textit{Nu}_m$ and the Rayleigh–Darcy number $ \textit{Ra}$ . Our experimental data clearly resolve the conductive regime and the onset of convection, with the value of onset Rayleigh–Darcy number agreeing well with the theoretical prediction of $ \textit{Ra}_c=4\pi ^2$ (Lapwood Reference Lapwood1948). In the moderate- $ \textit{Ra}$ range ( $4\pi ^2\leqslant Ra\lt 1000$ ), our results exhibit remarkable consistency with the DOB data, indicating that the experimental system can be effectively described by the DOB equations. Whereas for $ \textit{Ra}\gt 1000$ , our data undergo a gradual transition to the classical RBC scaling of $ \textit{Nu}_m\sim Ra^{0.3}$ , deviating significantly from those obtained by DOB simulations, which are close to $ \textit{Nu}_m\sim Ra$ . This deviation signifies the breakdown of the Darcy law approximation in the experiments. With increasing $ \textit{Ra}$ , the inertial effect – neglected in the DOB equations – becomes more and more important, eventually driving the system toward the classical RBC scaling. The transition from an RDC-like regime to an RBC-like regime also implies a drastic change in the flow structures and the dominant mechanism of energy dissipation.

In addition to Darcy’s law, another key assumption of DOB simulation is local thermal equilibrium (LTE) between the fluid and solid phases. In our experiments, the validity of LTE also depends on the Rayleigh–Darcy number: at small $ \textit{Ra}$ , the flow is weak, and the heat exchange between the fluid and solid phases is sufficiently fast to satisfy LTE. As $ \textit{Ra}$ increases, the flow becomes more intense, and the thermal relaxation between the two phases becomes less effective, so that LTE is expected to gradually break down at sufficiently high $ \textit{Ra}$ . The progressive departure from LTE may also contribute to the observed discrepancy between our experimental results and existing DOB simulations at high $ \textit{Ra}$ .

Figure 3.

(a) The porous medium Nusselt number $ \textit{Nu}_m$ as a function of Rayleigh–Darcy number $ \textit{Ra}$ . The solid lines show fittings with (3.1). The dashed black line and blue line show the fitting functions adopted from references Keene & Goldstein (Reference Keene and Goldstein2015) and Pirozzoli et al. (Reference Pirozzoli, De Paoli, Zonta and Soldati2021), respectively. The blue solid diamonds represent values of $ \textit{Ra}$ at which the centre-temperature measurements are conducted. (b) The local scaling exponent $\beta$ of the $ \textit{Nu}_m$ – $ \textit{Ra}$ curve as a function of $ \textit{Ra}$ . The solid lines correspond to those in panel (a). The data obtained from classical RBC are plotted with solid squares for benchmarking, where the $ \textit{Ra}_{\kern-1pt f}$ axis is located at the top of the figure.

3.2. Regimes of heat transport

For a quantitative analysis of the transitional behaviour of heat transport, we adopt the following $n$ th-order transition function (Zhang & Xia Reference Zhang and Xia2023b )

(3.1) \begin{align} Nu_m & = A\boldsymbol{\cdot }Ra^{\beta _n}\boldsymbol{\cdot }\left [1+\left (\frac {Ra_{t,1}}{Ra}\right )^4\right ]^{-\frac {\beta _1-\beta _2}{4}} {\cdots} \nonumber \\ & \quad \left [1+\left (\frac {Ra_{t,n-1}}{Ra}\right )^4\right ]^{-\frac {\beta _{n-1}-\beta _{n}}{4}}, \qquad n\geq 2, \end{align}

where $A$ , $\beta _{1\ldots n}$ , $ \textit{Ra}_{t,1 \ldots n-1}$ are $2n$ free fitting parameters. This function exhibits the following scalings: $ \textit{Nu}_m\sim Ra^{\beta _1}$ for $ \textit{Ra}\ll Ra_{t,1}$ ; $ \textit{Nu}_m\sim Ra^{\beta _i}$ for intermediate $ \textit{Ra}$ regimes ( $ \textit{Ra}_{t,i-1}\ll Ra\ll Ra_{t,i}$ ); and $ \textit{Nu}_m\sim Ra^{\beta _n}$ for $ \textit{Ra}\gg Ra_{t,n-1}$ . Thus, it allows $ \textit{Nu}_m$ to follow different power-law scalings (with exponent $\beta _1$ to $\beta _n$ ) across successive $ \textit{Ra}$ regimes.

We fitted our experimental data over the range $4\pi ^2\lt Ra\leqslant 2.62\times 10^5$ using (3.1) with $n$ up to 5. For comparison, the three-dimensional DOB data of Hewitt, Neufeld & Lister (Reference Hewitt, Neufeld and Lister2014) (over $4\pi ^2\lt Ra\lt 1000$ ) were fitted with a third-order form of (3.1). The fitting curves are plotted as solid lines in figure 3(a) and the corresponding fitting parameters are listed in table 2.

Table 2.

Fitting parameters of the transition functions.

Table 2. Long description

A table comparing fitting parameters of transition functions. The table has two rows and nine columns. Column headers are A, Ra_t_1, Ra_t_2, Ra_t_3, Ra_t_4, beta_1, beta_2, beta_3, beta_4, and beta_5. Row labels are Present study and 3-D DOB. Row 1: A, 3.59; Ra_t_1, 87.7; Ra_t_2, 455; Ra_t_3, 1.05 x 10^3; Ra_t_4, 1.05 x 10^4; beta_1, 1.29; beta_2, 0.491; beta_3, 1.522; beta_4, 0.589; beta_5, 0.335. Row 2: A, 0.01; Ra_t_1, 89.3; Ra_t_2, 351; Ra_t_3, ---; Ra_t_4, ---; beta_1, 1.17; beta_2, 0.313; beta_3, 1.03; beta_4, ---; beta_5, ---.

As shown in figure 3(a), the fifth-order transition function provides a precise quantitative description of the present experimental data across the entire $ \textit{Ra}$ range, validating its suitability for describing the multi-stage transitional heat transport in PMC. For a clearer view of the transition, we plot in figure 3(b) the corresponding local scaling exponent $\beta$ . The solid curves are extracted from the derivatives of the fitted equations, while the diamonds are determined using a power-law fitting to a sliding window of three adjacent experimental data points. The excellent agreement between these two approaches justifies the rationality of selecting a fitting function with $n=5$ .

With the above quantitative description of $\beta$ – $ \textit{Ra}$ relationship, we can identify five distinct regimes, marked by different colours in figures 3(a) and 3(b). The critical Rayleigh–Darcy numbers demarcating these regimes are denoted as $ \textit{Ra}_{\textit{II}}$ , $ \textit{Ra}_{\textit{III}}$ , $ \textit{Ra}_{\textit{IV}}$ and $ \textit{Ra}_{{V}}$ . The first critical point $ \textit{Ra}_{\textit{II}}$ marks the onset of convection and is consistent with the theoretical prediction of $4\pi ^2$ . The second and third critical points, $ \textit{Ra}_{\textit{III}}$ and $ \textit{Ra}_{\textit{IV}}$ were determined from the local extrema of ${\rm d}\beta /d(\log Ra)$ . To determine $ \textit{Ra}_{{V}}$ , we extended the tangent to $\beta (Ra)$ at the last minimum of ${\rm d}\beta /d(\log Ra)$ until it intersected the asymptotic value $\beta =\beta _5$ , and defined this intersection as $ \textit{Ra}_{{V}}$ . The details are provided in Appendix B. It is also worth mentioning that the values $ \textit{Ra}_{\textit{III}}$ , $ \textit{Ra}_{\textit{IV}}$ and $ \textit{Ra}_{{V}}$ are not identical to the fitting parameters of (3.1).

For $0\leqslant Ra\leqslant Ra_{\textit{II}}=4\pi ^2$ , heat transport is purely conductive, with $ \textit{Nu}_m=1$ and $\beta \approx 0$ . Convection sets in at $ \textit{Ra}_{\textit{II}}$ . The exponent $\beta$ jumps sharply to a value greater than unity and then gradually decreases to approximately 0.5 as $ \textit{Ra}$ increases to around $ \textit{Ra}_{\textit{III}}$ . As $ \textit{Ra}$ exceeds $ \textit{Ra}_{\textit{III}}$ , $\beta$ exhibits a significant jump to a value greater than one. Beyond $ \textit{Ra}_{\textit{IV}}$ , the local scaling exponent again shows a mild decrease. Finally, for $ \textit{Ra}$ greater than $ \textit{Ra}_{\text{V}}$ , $\beta$ approaches a value of approximately 0.335.

Among the identified regimes, regime I corresponds to the conductive state before the onset of convection. The characteristics of regime II and regime III are consistent with the results obtained in RDC, as evidenced by the close agreement between the experimental and DOB $\beta$ – $ \textit{Ra}$ curves in the corresponding $ \textit{Ra}$ range (figure 3 b). This indicates that the system’s behaviour in these two regimes is primarily governed by Darcy’s law, and also justifies them as the convection regime (II) and oscillation regime (III) (we will confirm this point in § 3.3). In contrast, the asymptotic scaling in regime V closely aligns with the heat transport scaling of the classical RBC system, indicating that the system’s dynamics is now dominated by inertial effects and viscous dissipation. Therefore, regime V is identified as the classical RBC regime. Lastly, regime IV represents the transitional process where the scaling exponent $\beta$ decreases systematically from RDC-like scaling toward the RBC-like asymptotic value. Physically, this regime corresponds to the system’s transition from a Darcy-drag-dominated to an inertia-dominated dynamics. We therefore identify it as the transition regime, which serves as a bridge connecting the RDC and RBC systems.

To summarise, the application of the transition function (3.1) has enabled a quantitative identification of five consecutive regimes: I. conduction regime; II. convection regime; III. oscillation regime; IV. transition regime; and V. classical RBC regime.

It is noteworthy that the confinement parameter $\varLambda$ , proposed by Noto, Letelier & Ulloa (Reference Noto, Letelier and Ulloa2024) and first introduced to PMC by Schwendener et al. (Reference Schwendener, Noir, Latt, Coreixas and Kong2026), provides an alternative perspective for interpreting the flow transition. Defining $\varLambda =\delta _{\textit{th}}/l_p=(H/l_p)/(2\textit{Nu}_m)$ as the ratio of the thermal boundary-layer thickness to the pore size, we find that the critical Rayleigh–Darcy number $ \textit{Ra}_{\textit{IV}}$ corresponds to $\varLambda \approx 1.5$ , i.e. $\varLambda =O(1)$ . This suggests that the transition from regime III to regime IV occurs when the thermal boundary layer becomes comparable to the pore scale, and can thus be interpreted as a consequence of plume-scale confinement relaxation.

3.3. Flow states at different regimes

Figure 4.

(a–f) The time series of the experimentally measured temperature fluctuation $T_c-\langle T_c\rangle$ at the cell centre for different values of $ \textit{Ra}$ . The horizontal axis is normalised by the free-fall time $t_{ff}=\sqrt {H/(\alpha g\Delta T)}$ . (g–l) The PSDs of the corresponding time series in (a–f). The black line denotes the $f^{-5/3}$ power law.

In this subsection, we present results from centre-temperature ( $T_c$ ) measurements, which shed light on the evolution of the flow state with increasing $ \textit{Ra}$ .

Figure 4(a–f) shows the time series of $T_c-\langle T_c\rangle$ measured at six representative values of $ \textit{Ra}$ – marked as blue solid diamonds in figure 3(a). At $ \textit{Ra}=261$ (within regime II), the centre-temperature time series (figure 4 a) remains nearly constant, with its corresponding power spectrum distribution (PSD) exhibiting white noise characteristics (figure 4 g). This time-independent local temperature behaviour confirms that the flow is in a steady convection state. As $ \textit{Ra}$ increases to $ \textit{Ra}=655$ , the system enters the oscillation regime (regime III). The centre temperature displays low-frequency, small-amplitude oscillations (figure 4 b), which is also reflected in its PSD (figure 4 h). A further increase of $ \textit{Ra}$ to $ \textit{Ra}=1.04\times 10^3$ (still within regime III) excites oscillations with multiple harmonic frequencies (see figures 4 c and 4 i).

When the system is in the transition regime ( $ \textit{Ra}=5.21\times 10^3$ ), the time series of centre temperature shows irregular, intermittent fluctuations (see figure 4 d). The corresponding PSD (figure 4 j) shows a smooth broadband distribution, indicating a transition toward turbulence. At $ \textit{Ra}=3.27\times 10^4$ , the time series of $T_c$ displays more intense fluctuations (figure 4 e), with its PSD (figure 4 k) exhibiting a power-law scaling close to $k^{-5/3}$ – a hallmark of turbulence (Sano, Wu & Libchaber Reference Sano, Wu and Libchaber1989; Niemela et al. Reference Niemela, Skrbek, Sreenivasan and Donnelly2000; Zhou & Xia Reference Zhou and Xia2001).

Notably, at $ \textit{Ra}=1.31\times 10^5$ , the intensity of the fluctuation of the centre temperature is significantly lower than that at $ \textit{Ra}=3.27\times 10^4$ . This is attributed to the fact that measurements for these two $ \textit{Ra}$ were taken in cells with different aspect ratios ( $\varGamma _{xy}$ ) (see table 1): the time series in figures 4(e) and 4(k) were obtained in a cell with $\varGamma _{xy}=1$ , while those of figures 4(f) and 4(l) were measured in a cell with $\varGamma _{xy}=1/2$ instead. In convective flows, domains with different aspect ratios usually develop distinct large-scale flow structures (Xi & Xia Reference Xi and Xia2008b ; van der Poel et al. Reference van der Poel, Stevens, Sugiyama and Lohse2012; Wang et al. Reference Wang, Verzicco, Lohse and Shishkina2020). Therefore, for these two sets of local temperature measurements, the ‘geometric centre’ may correspond to regions with different large-scale flow characteristics, leading to the seemingly singular temperature fluctuations behaviour. This physical picture is validated by preliminary results of flow visualisation, which will be presented later in § 3.5.

3.4. Effect of Darcy number $ \textit{Da}$

As mentioned in § 1, one of the central advantages of employing 3D-printed matrices as the porous medium is that it enables precise control of the Darcy number $ \textit{Da}$ . In this subsection, we will explore in detail how the Darcy number affects the transitions between different flow regimes. To achieve this, we conducted a series of additional experiments with different combinations of $l_p$ and $H$ , covering a Darcy number range of $6.18\times 10^{-7}\leqslant Da \leqslant 1.21\times 10^{-5}$ , while maintaining a constant porosity of $\varphi = 0.5$ .

Figure 5(a) presents the compensated medium Nusselt number, $ \textit{Nu}_m/Ra$ , as a function of $ \textit{Ra}$ for different $ \textit{Da}$ values. In general, the five datasets exhibit nearly identical shapes but shift toward the lower right as the Darcy number decreases. For quantitative analysis, each dataset (within a Rayleigh–Darcy number range of $100\lt Ra\lt 2\times 10^4$ ) was fitted using the third-order transition function ((3.1), $n=3$ ). To mitigate the adverse effects of overfitting, we first fitted the dataset of $Da=1.21\times10^{-5}$ to obtain $\beta _1=0.417$ . This value of $\beta _1$ was subsequently fixed for all other datasets during fitting, while the rest of the parameters remained unconstrained. The fitting results are plotted as solid curves in figure 5(a), and the fitting parameters are summarised in table 3. The two critical Rayleigh–Darcy numbers, $ \textit{Ra}_{\textit{III}}$ and $ \textit{Ra}_{\textit{IV}}$ (marked by solid circles and squares), were determined for each $ \textit{Da}$ using the same method detailed in § 3.2. Specifically, we extracted the local scaling exponent $\beta$ as a function of $ \textit{Ra}$ from the derivative of the third-order transition function and identified the critical $ \textit{Ra}$ values by finding the extrema of ${\rm d}\beta /d(\log Ra)$ . This approach enables a direct comparison of the transitional behaviour across different $ \textit{Da}$ .

Table 3.

Fitting parameters of the transition functions for different Darcy numbers.

Table 3. Long description

The table presents fitting parameters of transition functions for various Darcy numbers. It contains six rows and eight columns. The columns are labeled as Da, A, Ra_I,1, Ra_I,2, β1, β2, β3, Ra_III, and Ra_IV. The rows list different Darcy numbers and their corresponding fitting parameters. Row 1: Da, 1.21 times 10 to the power of -5; A, 0.353; Ra_I,1, 5.45 times 10 to the power of 2; Ra_I,2, 8.65 times 10 to the power of 2; β1, 0.417; β2, 1.96; β3, 0.59; Ra_III, 5.50 times 10 to the power of 2; Ra_IV, 8.60 times 10 to the power of 2. Row 2: Da, 5.39 times 10 to the power of -6; A, 0.553; Ra_I,1, 1.04 times 10 to the power of 3; Ra_I,2, 2.23 times 10 to the power of 3; β1, 0.417; β2, 1.47; β3, 0.55; Ra_III, 1.05 times 10 to the power of 3; Ra_IV, 2.30 times 10 to the power of 3. Row 3: Da, 2.47 times 10 to the power of -6; A, 0.318; Ra_I,1, 1.90 times 10 to the power of 3; Ra_I,2, 2.09 times 10 to the power of 3; β1, 0.417; β2, 7.04; β3, 0.63; Ra_III, 1.91 times 10 to the power of 3; Ra_IV, 2.10 times 10 to the power of 3. Row 4: Da, 1.35 times 10 to the power of -6; A, 0.501; Ra_I,1, 2.49 times 10 to the power of 3; Ra_I,2, 3.17 times 10 to the power of 3; β1, 0.417; β2, 3.29; β3, 0.59; Ra_III, 2.50 times 10 to the power of 3; Ra_IV, 3.10 times 10 to the power of 3. Row 5: Da, 6.18 times 10 to the power of -7; A, 0.376; Ra_I,1, 3.00 times 10 to the power of 3; Ra_I,2, 5.95 times 10 to the power of 3; β1, 0.417; β2, 1.50; β3, 0.62; Ra_III, 3.02 times 10 to the power of 3; Ra_IV, 6.05 times 10 to the power of 3.

Figure 5.

(a) Compensated plot of the experimentally measured media Nusselt number $ \textit{Nu}_m/Ra$ versus $ \textit{Ra}$ for different values of $ \textit{Da}$ . The solid lines represent fitting curves of (3.1) with $n=3$ . (b) The phase diagram of different heat transport regimes in the $ \textit{Ra}$ – $ \textit{Da}$ space. The solid symbols in both panel (a) and panel (b) present the critical Rayleigh–Darcy numbers.

For the Darcy number range explored in this study, it is seen that the critical Rayleigh–Darcy numbers ( $ \textit{Ra}_{\textit{III}}$ and $ \textit{Ra}_{\textit{IV}}$ ) – marking the transitions from regime II to regime III and from regime III to regime IV – both shift to higher values with decreasing $ \textit{Da}$ (i.e. reduced permeability). Furthermore, both $ \textit{Ra}_{\textit{III}}$ (onset of time-dependent flow) and $ \textit{Ra}_{\textit{IV}}$ (onset of the transition to RBC) exhibit power-law dependencies on $ \textit{Da}$ with $ \textit{Ra}_{\textit{III}}\sim Da^{-0.584}$ and $ \textit{Ra}_{\textit{IV}}\sim Da^{-0.586}$ , respectively (see the black solid lines in figure 5 b). The exponents are nearly identical, implying a consistent mechanism governing both transitions.

Previous studies of pore-scale DNS suggest that the transition from regime III to regime IV occurs when the thermal boundary-layer thickness $\delta _{\textit{th}}$ becomes comparable to the pore size $l_p$ (Liu et al. Reference Liu, Jiang, Chong, Zhu, Wan, Verzicco, Stevens, Lohse and Sun2020), which corresponds to a confinement parameter of order unity $\varLambda = O(1)$ . Building on this hypothesis, if we make a first-order assumption that the $ \textit{Nu}_m$ – $ \textit{Ra}$ relationship in both regime II and regime III follows a single power-law scaling $ \textit{Nu}_m\sim Ra^{\bar \beta }$ , and consider the fact that $ \textit{Da}=K/H^2\sim l_p^2/H^2$ and $\delta _{\textit{th}}/H=1/2Nu_m$ , we can obtain the following relationship for the transitional Rayleigh–Darcy number associated with $\delta _{\textit{th}}\sim l_p$ :

(3.2) \begin{align} Ra_{\textit{IV}}\sim Da^{-1/(2\bar \beta )}. \end{align}

The experimentally measured scaling exponent $ \textit{Ra}_{\textit{IV}}\sim Da^{-0.586}$ yields $\bar \beta = 0.853$ , which is consistent with the measured local scaling exponents (figure 3 b). This result further validates the hypothesis that the transition from regime III to regime IV is associated with the cross-over between the thermal boundary-layer thickness and the pore size.

In addition, the critical Rayleigh–Darcy number $ \textit{Ra}_{{V}}$ for $ \textit{Da}=5.39\times 10^{-6}$ was also obtained from its $ \textit{Nu}_m(Ra)$ fitting curve. The boundary between regimes IV and V in the phase diagram (dashed line in figure 5 b) was then estimated by connecting the two solid diamonds with a power law. This power-law dependence of $ \textit{Ra}_{{V}}$ on $ \textit{Da}$ , $ \textit{Ra}_{{V}}\sim Da^{-0.220}$ , differs significantly from those of $ \textit{Ra}_{\textit{III}}$ and $ \textit{Ra}_{\textit{IV}}$ , indicating a distinct physical mechanism governing the transition to the classical RBC regime.

With the above, we can now construct a phase diagram of the flow states in the $ \textit{Ra}$ – $ \textit{Da}$ parameter space, as shown in figure 5(b). This diagram clarifies the role of the Darcy number: it shows a significant influence on the transition between regimes II, III and IV. At low $ \textit{Da}$ (dense porous media), the system exhibits characteristic RDC behaviour (regimes II and III) over a relatively wide range of $ \textit{Ra}$ ; in contrast, at high $ \textit{Da}$ (sparse porous media), the system evolves into the transition regime (IV) and the classical RBC regime (V) at a substantially lower $ \textit{Ra}$ . Although the critical point that separates regime III and regime IV can be rationalised by equalising the thermal boundary-layer thickness $\delta _{\textit{th}}$ and the pore size $l_p$ , the physical mechanisms governing the transitions among other regimes remain unresolved and merit future investigation. Nevertheless, our results successfully reconcile the seemingly contradictory findings from DOB simulations, packed-bead experiments and pore-scale DNS.

Beyond the five flow regimes discussed above, there is an alternative classification of flow regimes based on the relative importance of the inertial effects: the Darcy regime and the Forchheimer regime (Nield & Bejan Reference Nield and Bejan2017). In the former, flow is governed by Darcy’s law, with the buoyancy force balanced by the Darcy drag, resulting in the linear heat transport scaling, $ \textit{Nu}_m\sim Ra$ . As flow intensity increases, the linearity of Darcy’s law gradually breaks down, and a modification with Forchheimer drag becomes necessary (Joseph, Nield & Papanicolaou Reference Joseph, Nield and Papanicolaou1982). In the Forchheimer regime, a balance between buoyancy and Forchheimer drag yields $ \textit{Nu}_m\sim Ra^{1/2} \textit{Pr}_p^{1/2}$ , where $ \textit{Pr}_p= \textit{Pr}_fc_F^{-1}Da^{-1/2}\lambda$ is the modified Prandtl number with $c_F$ being the Forchheimer coefficient (Wang & Bejan Reference Wang and Bejan1987). The transition point from the Darcy regime to the Forchheimer regime can be estimated by equating the two scaling relationships above, yielding $ \textit{Ra}_{D-F}\sim \textit{Pr}_p$ . Substituting the definition of the modified Prandtl number, we can obtain $ \textit{Ra}_{D-F}\sim Da^{-1/2}$ , which is generally consistent with the transition from regime III to regime IV identified in this study.

3.5. The flow structures in regimes IV and V

Figure 6.

(a–d) Experimentally measured velocity fields at the centre $xy$ -plane ( $z=1/2W$ ) for (a) $ \textit{Ra}=2.12\times 10^3$ , (b) $ \textit{Ra}=1.06\times 10^4$ and (c) $ \textit{Ra}=3.35\times 10^4$ in a $\varGamma _{xy}=1$ cell; (d) $ \textit{Ra}=1.34\times 10^5$ in a $\varGamma _{xy}=1/2$ cell. (e) Pore-scale Reynolds number $ \textit{Re}_p$ based on maximum or volume-averaged root-mean-square velocity as a function of Rayleigh–Darcy number $ \textit{Ra}$ . The solid lines represent the power-law fit to the experimental data as a first-order approximation, and the dashed lines show the extrapolation of the fitted function. Error bars at both ends of the fitted curve represent a 90 % confidence interval. The Darcy number is fixed at $ \textit{Da}=1.21\times 10^{-5}$ .

Figure 6. Long description

The image contains four vector field plots and one log-log plot. Panel A, B, and C show experimentally measured velocity fields at the center-plane for different Rayleigh-Darcy numbers in a cell. Each panel displays a vector field plot with x and y axes labeled in millimeters, and velocity vectors color-coded according to their magnitude in millimeters per second. Panel D shows a similar vector field plot but for a different cell configuration. Panel E presents a log-log plot of the pore-scale Reynolds number based on maximum or volume-averaged root-mean-square velocity as a function of the Rayleigh-Darcy number. The x-axis represents the Rayleigh-Darcy number, and the y-axis represents the Reynolds number. Solid lines represent power-law fits to the experimental data, while dashed lines show the extrapolation of the fitted function. Error bars at both ends of the fitted curve represent a 90 percent confidence interval.

As mentioned above, another key advantage of the present experimental approach is that it enables visualisation of the flow field within the porous media. In this subsection, we present preliminary results of the global flow structure measured via PIV. The measurements were conducted in the vertical $xy$ -plane located at $z=1/2W$ of cell PMC2 and cell PMC3 (see table 1). The investigated values of the Rayleigh–Darcy number $ \textit{Ra}=2.12\times 10^3$ , $1.06\times 10^4$ , $3.35\times 10^4$ and $1.34\times 10^5$ cover the transition regime and classical RBC regime.

Figure 6(a–d) shows the measured velocity fields (averaged over about one minute) within the porous media. The flow is seen to evolve from several horizontally stacked convection rolls to a single large-scale circulation. The number of rolls $n$ , with a typical horizontal length scale of $l=L/n$ , appears to decrease with increasing $ \textit{Ra}$ . On the other hand, the non-dimensional horizontal length scale of the flow structures – $l/H=(L/n)/H=\varGamma _{xy}/n$ – exhibits a monotonic non-decreasing trend – it increases from approximately 1/3 (figure 6 a) to 1/2 (figure 6 b), and then stays constant as $ \textit{Ra}$ further increases (see figure 6 b–d). This is in contrast to DOB simulations, which report a negative power-law scaling between $l/H$ and $ \textit{Ra}$ (Wen, Corson & Chini Reference Wen, Corson and Chini2015; De Paoli et al. Reference De Paoli, Pirozzoli, Zonta and Soldati2022; Zhu et al. Reference Zhu, Fu and De Paoli2024). This discrepancy originates from the fact that DOB simulations rely on the Darcy law (neglecting inertia), whereas the present PIV measurements are conducted in regime IV and regime V (see figure 6 e), where the Darcy law is no longer valid. Notably, even in the classical RBC regime, the global flow structures obtained by the present study (see figures 6 c and 6 d) are apparently different from those measured in a classical RBC system, where a single large-scale circulation (for cells with unity aspect ratio) (Xia, Sun & Zhou Reference Xia, Sun and Zhou2003; Xi & Xia Reference Xi and Xia2008a ) or two vertically stacked convection rolls (for cell with aspect ratio one half) (Stringano & Verzicco Reference Stringano and Verzicco2006; Xi & Xia Reference Xi and Xia2008b ) are reported.

Furthermore, we calculated the pore-scale Reynolds number $ \textit{Re}_p$ based on the volume-averaged root-mean-square (r.m.s.) velocity, $ \textit{Re}_{p,\langle u_{\textit{rms}}\rangle _V}=\langle u_{\textit{rms}}\rangle _Vl_p/\nu_{\kern-1pt f}$ , and that based on the maximum r.m.s. velocity $ \textit{Re}_{p,u_{\textit{rms},\textit{max}}}= u_{\textit{rms},\textit{max}}l_p/\nu_{\kern-1pt f}$ . Figure 6(e) shows $ \textit{Re}_p$ as a function of the Rayleigh–Darcy number $ \textit{Ra}$ . For all explored $ \textit{Ra}$ , both definitions of pore-scale Reynolds number are considerably larger than unity, indicating that Darcy’s law is invalid – consistent with our heat transport flow regime classification. In addition, both $ \textit{Re}_{p,\langle u_{\textit{rms}}\rangle _V}$ and $ \textit{Re}_{p,u_{\textit{rms},\textit{max}}}$ exhibit power-law relationships with $ \textit{Ra}$ , with best fits yielding $ \textit{Re}_{p,u_{\textit{rms},\textit{max}}}\sim Ra^{0.648}$ and $ \textit{Re}_{p,\langle u_{\textit{rms}}\rangle _V}\sim Ra^{0.678}$ , respectively. It is worth mentioning that $ \textit{Re}_p$ is evaluated based on PIV measurements within the porous structure, where the optical interrogation window is limited by the solid matrix, leading to considerable uncertainties. Therefore, the single power-law fits for the $ \textit{Re}_p$ – $ \textit{Ra}$ relationships only serve as a first-order approximation. Nevertheless, the obtained scaling exponents are slightly higher than the 1/2 value predicted for the Forchheimer regime. And the intersection points of the two extrapolated curves with $ \textit{Re}_p=1$ are close to the onset of convection and oscillation, respectively.

It is also worth mentioning that mechanical dispersion (De Paoli Reference De Paoli2023), an inherent effect in PMC, may also influence the transitions between different flow regimes. Mechanical dispersion arises from spatial variability in pore-scale velocities and is responsible for the broadening of thermal plumes induced by inter-pore and intra-pore advection. A quantitative characterisation of mechanical dispersion would require independent measurement of the mechanical dispersion coefficient, or introducing an additional dimensionless parameter – the dispersion Rayleigh–Darcy number $ \textit{Ra}_d$ (Liang et al. Reference Liang, Wen, Hesse and DiCarlo2018). However, this demands high-resolution flow field and temperature field measurements, and therefore is beyond the scope of the present experimental study.