Pore-scale statistics of temperature and thermal energy dissipation rate in turbulent porous convection

In the pore-scale method, individual pore geometry is directly resolved, thus, the governing equations for fluid flows and heat transfer in the pores are the Navier-Stokes equations with Boussinesq approximation:

Here, the subscript $f$ denotes the fluid phase. $\mathbf{u}_{f}$, $P$, and $T_{f}$ are the fluid velocity, pressure, and temperature in the pores, respectively. $\rho_{0}$ and $T_{0}$ are reference density and temperature, respectively. $g$ is the gravity value and $\hat{\mathbf{y}}$ is the unit vector in the vertical direction. In Eq. (1c), all the transport coefficients (i.e., $\nu_{f}$, $\alpha_{f}$, $\beta$) are assumed to be constants. Using the nondimensional group

Equation (1c) can be rewritten in dimensionless form as

In the following, for convenience, we will drop the superscript star ($*$) to denote a dimensionless variable. The dimensionless parameters of the Rayleigh number ($Ra$) and the Prandtl number ($Pr$) are defined as

We adopt the lattice Boltzmann (LB) method [33, 34, 35] as the numerical tool for direct numerical simulation of turbulent thermal convection in the pores. In the LB method, to solve Eqs. (1ca) and (1cb), the evolution equation of density distribution function is written as

To solve Eq. (1cc), the evolution equation of temperature distribution function is written as

Here, $f_{i}$ and $g_{i}$ are the density and temperature distribution functions, respectively. $\mathbf{x}$ is the fluid parcel position, $t$ is the time, $\delta_{t}$ is the time step. $e_{i}$ is the discrete velocity along the $i$ th direction. $\mathbf{M}$ is a $9\times 9$ orthogonal transformation matrix based on the D2Q9 discrete velocity model; $\mathbf{N}$ is a $5\times 5$ orthogonal transformation matrix based on the D2Q5 discrete velocity model. The equilibrium moments $\mathbf{m}^{(\text{eq})}$ in Eq. (5) are

The equilibrium moments $\mathbf{n}^{(\text{eq})}$ in Eq. (6) are

where $a_{T}$ is a constant determined by the thermal diffusivity as $a_{T}=20\sqrt{3}\kappa_{f}-4$. The relaxation matrix $\mathbf{S}$ is $\mathbf{S}=\operatorname{diag}\left(s_{\rho},s_{e},s_{s},s_{j},s_{q},s_{j},s_{q},s_{\nu},s_{\nu}\right)$, and the kinematic viscosity of the fluids is calculated as $v_{f}=\left(s_{v}^{-1}-0.5\right)/3$. The relaxation matrix $\mathbf{Q}$ is given by $\mathbf{Q}=\operatorname{diag}\left(0,q_{k},q_{k},q_{e},q_{v}\right)$, where $q_{\kappa}=3-\sqrt{3}$ and $q_{e}=q_{\nu}=4\sqrt{3}-6$.

The macroscopic fluid variables of density $\rho_{f}$, velocity $\mathbf{u}_{f}$ and temperature $T_{f}$ are calculated as $\rho_{f}=\sum_{i=0}^{8}f_{i}$, $\mathbf{u}_{f}=\left(\sum_{i=0}^{8}\mathbf{e}_{i}f_{i}+\mathbf{F}/2\right)/\rho_{f}$ and $T_{f}=\sum_{i=0}^{4}g_{i}$, respectively. More numerical details on the LB method and validation of the in-house code can be found in our previous work [36, 37, 38].

We assume the solid matrix is impermeable to both fluid and heat flux. At the fluid-solid interface, the no-slip velocity boundary conditions can be described as $\mathbf{u}_{f}=0$; while the adiabatic temperature boundary conditions can be described as $\partial_{\mathbf{n}}T_{f}=0$. In the LB method, the-above fluid-solid interface conditions can be mimic by the bounce-back rules for the density and temperature distribution functions as $f_{\bar{i}}\left(\mathbf{x}_{f},t+\delta_{t}\right)=f_{i}^{*}\left(\mathbf{x}_{f},t\right)$ and $g_{\bar{i}}\left(\mathbf{x}_{f},t+\delta_{t}\right)=g_{i}^{*}\left(\mathbf{x}_{f},t\right)$, respectively.

We artificially construct the porous structure via randomly placing square cylinders of length $d$ in the RB convection cell, as illustrated in Fig. 1(a). In Fig. 1(b), we further present an enlarged view of the porous region as that in Fig. 1(a), and we illustrate the directional average method [39, 40] to calculate pore size distribution (PSD). At each fluid point, we start counting the pore length along with specified directions until reaching a solid point; then, the pore diameter is obtained by averaging the pore length in all given eight directions. It should be noted that the directional average method provides an approximate assessment of pore size, and it comes with inherent limitations. This method determines the average span of the pore spaces from a certain point in all possible directions (i.e., eight discrete directions in this work) until an obstacle or the boundary of the domain is encountered. When assessing fluid points near the boundary, we exclude directions beyond the domain boundary. The absence of obstacles in certain directions could lead to an unusually large distance, causing an overestimation of pore size. Nevertheless, the method retains its applicability in this work, as the issues related to boundary effects or statistical anomalies can be mitigated by employing a large computational domain. In addition, it is particularly beneficial when our goal is to determine the average pore sizes across the porous medium. Figure 1(c) shows the probability density functions (PDFs) of calculated pore size for five different realizations of the porous structure at the same porosity of $\phi=0.86$. Here, the pore size is normalized by the cylinder length $d$. We can see that despite different realizations of the porous structure generation, the pore size distributions are roughly the same at the same porosity. In Fig. 1(d), we compare the pore size for porous structure with different porosities of $\phi=0.86$, 0.90 and 0.94. Overall, the pore size exhibits unimodal distribution, and it increases with the increase of porosity.

To determine the permeability $K$ of the porous structure, we performed another set of pore-scale simulations when the flow is isothermal and in the Darcy regime. As illustrated in Fig. 2(a), fluid flows through porous media is driven by a small pressure difference $(P_{in}-P_{out})/P_{ref}=2\times 10^{-4}$ either in lateral or longitudinal direction, such that flow is sufficient slow. Following the Darcy’s law [41], we calculate the permeability tensor as [12]

Here, the intrinsic phase average is defined as $\langle\psi_{f}\rangle^{f}=(1/V_{f})\int_{V_{f}}\psi_{f}dv$, and the superficial phase average is defined as $\langle\psi_{f}\rangle=(1/V)\int_{V_{f}}\psi_{f}dv$. $V_{f}$ denotes the volume of the fluid phase within the representative volume $V$, and $\psi_{f}$ is a quantity associated with the fluid phase. We also check the pore-scale Reynolds number $Re_{D}=\langle u_{f}\rangle^{f}D/\nu$ satisfies $Re_{D}\leq O(1)$, thus, the condition to apply Darcy’s law is guaranteed. In Fig. 2(b), we show the permeabilities in lateral and longitude directions (i.e., $K_{xx}$ and $K_{yy}$), respectively, which are further presented in terms of dimensionless Darcy number $Da=K/L^{2}$. Here, the characteristic length $L$ is chosen as convection cell size. We can see that the lateral and longitudinal permeabilities are generally the same, suggesting the artificially constructed porous structure is homogeneous. Besides, we compare the calculated permeability with empirical Blake-Kozeny-Carman relation [14]

We can see from Fig. 2(b), the empirical relation overestimates the permeability for the sparse porous media, and the deviations increase with increasing the porosity. For the investigated porosity range $0.86\leq\phi\leq 0.98$, the corresponding $Da$ range $10^{-4}<Da<10^{-2}$.

We consider a two-dimensional (2D) porous RB convection cell with size $L=H$. The top and bottom walls of the cell are kept at a constant cold and hot temperature, respectively; the other two vertical walls are adiabatic. All four walls impose no-slip velocity boundary conditions. We provide simulation results at Prandtl number of $Pr=5.3$ and $0.7$ (i.e., corresponding to the working fluid of water and air, respectively) and a fixed $Ra=10^{9}$. The porosity ($\phi$) range $0.86\leq\phi\leq 0.98$, and the corresponding porous Rayleigh number ($Ra^{*}=Ra\cdot Da$) range $10^{5}<Ra^{*}<10^{7}$, suggesting vigorous convection in porous media [16]. In addition, we show the scaling of the global quantities on one of the control parameters $Ra$ (for $10^{6}\leq Ra\leq 10^{9}$), while the porosity is fixed as $\phi=0.86$ and Prandtl number is fixed as $Pr=5.3$ and $0.7$. A total of 100 simulations were carried out for porous convection with impermeable solid matrix, and tabulated values on the results are listed in the Appendix. For the canonical RB convection, the mesh size of the convection cell is $1024$ l.u.$\times 1024$ l.u.; while for the porous RB convection, the mesh size is even finer with $1200$ l.u.$\times 1200$ l.u. The resolution for the cylinder length $d$ is 40 l.u., and the minima gap between two cylinders is 40 l.u. Here, l.u. denotes the lattice length unit in the LB simulation [42].

For canonical RB convection of pure fluid, we verify the grid spacing $\Delta_{g}$ and time interval $\Delta_{t}$ is properly resolved by comparing with the Kolmogorov and Batchelor scales. Specifically, the Kolmogorov length scale [43] is estimated by the global criterion $\eta_{K}=\left[\nu^{3}/(\varepsilon_{u})_{V,t}\right]^{1/4}=HPr^{1/2}/[Ra(Nu-1)]^{1/4}$ , the Batchelor length scale [44] is estimated by $\eta_{B}=\eta_{K}Pr^{-1/2}$, and the Kolmogorov time scale [43] is estimated as $\tau_{\eta}=\sqrt{\nu/\left\langle\varepsilon_{u}\right\rangle_{V,t}}=t_{f}\sqrt{Pr/(Nu-1)}$. Here, $\varepsilon_{u}$ denotes the kinetic energy dissipation rates, and its global average can be related to the Nusselt number via the exact relation [19] $\left\langle\varepsilon_{u}\right\rangle_{V,t}=\nu^{3}Ra(Nu-1)/(H^{4}Pr^{2})$. Simulation results have shown that grid spacing satisfies the criterion of max $\left(\Delta_{g}/\eta_{K},\ \Delta_{g}/\eta_{B}\right)\leq 0.55$, which ensures spatial resolution; the time intervals are $\Delta_{t}\leq 0.00047\tau_{\eta}$, thus adequate temporal resolution is guaranteed. In addition, we validate our results by comparing the Nusselt and Reynolds number with those obtained using the NEK5000 solver (version v19.0) [45], as well as previous results reported by Zhang et al. [25]. The tabulated values are presented in Table 1. Note that the small deviations in the average statistical value may be attributed to turbulence fluctuations. For porous RB convection, both the cylinder length ($d=40$ l.u.) and minima gap between the cylinders (i.e., $40$ l.u.) are much large than that of the boundary layer thickness (around $10$ l.u.), thus the pore space is adequately resolved.

A typical snapshot of an instantaneous temperature field in both a porous and a canonical RB convection cell is shown in Fig. 3, and the corresponding video can be viewed in the Supplemental Material [46]. We can see that the flow structure in the porous RB convection exhibits different patterns from that in the canonical RB convection. In the canonical RB convection [see Figs. 3(b) and 3(d)], the rising and falling thermal plumes self-organize into a well-defined large-scale circulation (LSC) that spans the size of the convection cell [47], and there exist counterrotating corner rolls. The temperature field is efficiently mixed in the convection cell, with the bulk temperature being almost a constant value of $(T_{hot}+T_{cold})/2$. In contrast, in the porous RB convection [see Figs. 3(a) and 3(c)], the flow structure is less coherent and the large-scale flow circulation is suppressed. The rising and falling plumes penetrate through the pore throat, resulting in less mixing of the temperature field. This flow pattern is similar to that observed in the study by Liu et al. [27], where the porous matrix was permeable to heat flux. However, in the current work, we assume that the solid porous matrix is impermeable to heat flux. When the solid porous matrix does not conduct heat, there is no thermal exchange between the solid and fluid phases, and the fluid phase’s ability to effectively interact with the solid phase is diminished compared to scenarios involving a thermally conductive solid porous matrix. Consequently, the plume dynamics are less coherent than in the previous study. In the Supplemental Material [46], we provide corresponding videos, which allow for a more detailed examination of the flow patterns and temperature field in the two types of convection cells.

To validate the above conjecture, we calculate the cross-correlation coefficient between vertical velocity $v$ and temperature $T$ along the mid-plane of the cell, given by

where $\sigma_{v,T}$ denotes the standard deviation of $v$ and $T$. We conducted simulations for two cases: one with solid porous matrix being permeable to heat flux (referred to as permeable heat flux, following the simulation settings reported by Liu et al. [27]), and the other is the solid porous matrix being impermeable to heat flux (referred to as impermeable heat flux). Five different realizations of the porous structure were considered at the same porosity $\phi=0.86$. Figures 4(a) and 4(b) shows the cross-correlation coefficient $R_{v,T}$ for both cases. We can observe that $R_{v,T}$ is generally lower for the impermeable heat flux case, which implies that the thermal plumes are less coherent. We further calculate the joint probability density distribution of the vertical velocity $v$ and temperature fluctuation $\delta T=T-(T_{hot}+T_{cold})/2$. In comparison to a thermally conductive solid porous medium, an impermeable medium reduces the correlation between vertical velocity and temperature. This effect is due to the inherent nonthermal conductivity property of the impermeable medium. When the solid porous matrix does not conduct heat, the fluid phase’s ability to interact effectively with the solid phase is diminished compared to scenarios with a heat-conductive solid porous matrix.

In Fig. 5, we show the time-averaged flow field (temperature field and streamlines), where we can see that the presence of the solid porous matrix disrupts the LSC (i.e., the tilted elliptical main roll at $Pr=5.3$ or the circular main roll at $Pr=0.7$ in the canonical RB cell), and the LSC shape becomes more irregular in the porous RB cell. We note the existence of solid surfaces with a no-slip boundary condition significantly affects the flow patterns, resulting in the differences in the smoothness of the streamlines between the porous RB cell and the canonical RB cell. Meanwhile, to address concerns regarding simulation convergence, we have checked temperature fields and streamlines averaged over different time intervals to demonstrate convergence (not shown here for clarity). The corner rolls in the porous cell are also suppressed due to the existence of the porous matrix. Overall, we expect such disruption would lead to much more complex substructures inside the LSC at some porosities. Specifically, when the porosity is too large, the solid porous matrix occupies only a small volume fraction in the convection cell and it will have minor effects on the flow structure; when the porosity is too small, detached plumes from the top and bottom walls will only penetrate through the porous throat, and the dense porous structure may prohibit the formation of the LSC. Once the substructures emerge inside the LSC, they contribute to an increased instability within the LSC, potentially inducing flow reversals in turbulent thermal convection [48]. Our simulations have confirmed this conjecture, as we indeed observe flow reversal in the porous RB convection at $Ra=10^{9}$ and $Pr=0.7$ with various porosities. Previous study by Sugiyama et al. [49] suggests that flow reversal is absent in the canonical 2D RB convection at the same $Ra$ and $Pr$ (i.e., $Ra=10^{9}$ and $Pr=0.7$), thus the solid porous matrix has a profound impact on the flow dynamics, and understanding these effects is essential for predicting and controlling convection in porous media.

We measure the global heat transport by the volume-averaged Nusselt number ($Nu$) as $Nu=1+\sqrt{PrRa}\left\langle v^{*}T^{*}\right\rangle_{V,t}$, while the global strength of the convection is measured by the Reynolds number ($Re$) as $Re=\sqrt{Ra/Pr}\sqrt{\left\langle u^{*2}+v^{*2}\right\rangle_{V,t}}$. Here, $\langle\cdot\rangle_{V,t}$ denotes the superficial phase and time average, the asterisk superscript (*) denote the dimensionless variables. At each porosity, we calculate the $Nu$ and $Re$ based on results from five different realizations of the porous structure. From Fig. 6(a), we can see that $Nu$ increases monotonously with the decreasing of porosity over $0.86\leq\phi\leq 0.98$ at $Pr=5.3$, while $Nu$ increases first and then decreases with the decreasing of porosity at $Pr=0.7$. The enhanced heat transfer efficiency with slightly decreasing porosity is attributed to strongly correlated velocity and temperature fields. However, further decreasing the porosity increases the impedance from the porous solid matrix on heat transfer. We hypothesize the competition of these two factors results in an optimal porosity value when heat transfer efficiency is maximized. For $Pr=0.7$, we observed this hypothesized optimal value at porosity around $0.90$, while for $Pr=5.3$, the optimal porosity value may be smaller than those under investigations (i.e., $\phi<0.86$), thus we did not observe such optimal value in our simulations. It should be noted that to ensure adequate resolution of the pore spaces, we only consider sparse porous media in our simulations; in addition, considering the observations alongside the presence of error bars, we recognize the necessity for caution when attempting to draw definitive conclusions regarding the existence of an optimal porosity. This is particularly crucial when accounting for the potential influence of varying Prandtl numbers. As for the global flow strength, we can see from Fig. 6(b) that $Re$ decreases monotonously with the decreasing of the porosity, which can be understood as the introduction of the porous solid matrix in the convection cell enhances flow resistance.

We also show the scaling of the global quantities, such as $Nu$ and $Re$, on one of the control parameters $Ra$ (for $10^{6}\leq Ra\leq 10^{9}$), while the porosity is fixed as $\phi=0.86$. We also provide $Nu$ and $Re$ in the canonical RB convection. Previously, Zhang et al. provided tabulated values of $Nu$ and $Re$ versus $Ra$ at $Pr=5.3$ and 0.7. Our simulation results on the canonical RB convection are in good agreement with those reported by Zhang et al. [25]. The data shown in Fig. 7 indicate that in the porous convection, the increase of $Nu$ and $Re$ gradually approaches the power-law relations $Nu\propto Ra^{0.30}$ and $Re\propto Ra^{0.59}$, consistent with previous results reported in the canonical RB convection [50, 51, 25]. At fixed $\phi$, the scaling behavior of $Nu$ and $Re$ with $Ra$ slightly deviates from that of the canonical RB convection when $Ra$ is smaller, suggesting that heat transfer and momentum exchange are not solely governed by the boundary layer.

Figure 8 shows the probability density functions (PDFs) of normalized temperature $(T-\mu_{T})/\sigma_{T}$ measured at three different heights. In both porous and canonical RB cells, the PDFs of temperature are left-skewed near the top region (i.e., $y=0.75H$) due to dominated cold falling plumes, right-skewed near the bottom region (i.e., $y=0.25H$) due to dominated hot rising plumes, and symmetric at mid-height (i.e., $y=0.5H$) as a result of comparable falling cold plumes and rising hot plumes. Meanwhile, despite these similarities, there are notable differences between the PDF profiles of the two cells. Specifically, in the canonical RB cell, the peak of the temperature PDF profiles exhibits a stretched exponential behavior, and the tails show a Gaussian behavior; in contrast, in the porous RB cell, the temperature PDF profiles are narrowed down, and the stretched exponential peaks are absent, indicating that porous media suppresses extreme temperature events in the convection cell.

To highlight the damping effect that arises from the presence of a porous structure, which impedes both hot and cold thermal plumes, we analyze the PDFs of the fluid temperature in two regions: near the solid porous matrix and away from it. Specifically, we consider fluids nodes with distances less than 5 l.u. from the porous matrix as the near region, and fluid nodes with distances more than 5 l.u. as the away region, as illustrated in Fig. 9(c). In Figs. 9(a) and 9(b), we plot the PDFs of the temperature obtained over the whole cell and over time in the above two regions. We can see that the PDFs of temperature in both regions exhibit a symmetric peak, indicating a comparable occurrence of hot and cold plumes in those two regions. However, the PDFs of temperature near the porous matrix have narrower tails compared to that away from the porous matrix. This narrower tail implies a reduced degree of small-scale intermittency of the temperature fluctuations near the porous matrix.

We now provide the averaged vertical profile of statistics of the temperature field to quantitatively describe the temperature distributions. We first calculate the averaged vertical profile of temperature $\langle T\rangle_{\mathbf{x},t}$ and mean square fluctuations $\langle\theta^{2}\rangle_{\mathbf{x},t}$ for various porosities, as shown in Fig. 10. Here, the fluctuations $\theta(\mathbf{x},t)=T(\mathbf{x},t)-\langle T\rangle_{\mathbf{x},t}(y)$; the average $\langle\cdot\rangle_{\mathbf{x},t}$ is calculated over time $t$ and along the horizontal line $\mathbf{x}$ in the fluid phase. Note that the average over the fluid phase is referred to as the intrinsic phase average, in contrast to the superficial phase average, which would encompass the entire porous media domain. For the porous RB convection, the vertical profiles are further averaged over five different realizations of the porous structure at the same porosity. This fivefold averaging is conducted to mitigate the statistical errors arising from the random distribution of the porous medium. We can see from Figs. 10(a) and 10(c), away from the top and bottom walls, the averaged vertical temperature profiles are almost a constant value of $(T_{hot}+T_{cold})/2$ in both canonical and porous RB cells. This finding suggests that the temperature field in the bulk region of the fluid is insensitive to the presence of the impermeable solid matrix. In contrast, Figs. 10(b) and 10(d) reveal that the averaged vertical profiles of mean-square fluctuations of the temperature are sensitive to the porosity. With the decreasing of porosity, the flow structure becomes less coherent, leading to an increase in the fluctuations of temperature. From the inset of Figs. 10(b) and 10(d), we can also measure the thickness of thermal boundary layer $\delta_{T}$ as the location of the peak value in the profile [22, 52], which is close to $H/(2Nu)$. For the canonical RB convection (i.e., $\phi=1.00$), the temperature fluctuation profile diverges between $Pr$ of 5.3 and 0.7. As shown in Fig. 10(b), the profile at $Pr=5.3$ exhibits two peaks around $y=0.65$ and $y=0.35$, different from the pattern observed at $Pr=0.7$ [see Fig. 10(d)]. This variation could be attributed to differences in the corner rolls at different Prandtl numbers. Upon comparing with Fig. 5(b), it becomes apparent that the enhancement in turbulent fluctuations coincides with heights of corner rolls that are situated diagonally, thereby leading to the emergence of the peaks observed in the bulk region of temperature fluctuation profile at $Pr=5.3$. These differences highlight the complexity of the system’s temperature fluctuations and plume dynamics in response to the Prandtl number.

We further calculate the averaged vertical profile of higher-order moments of temperature and plot the skewness $S_{\theta}(y)=\langle\theta^{3}\rangle_{\mathbf{x},t}/\langle\theta^{2}\rangle_{\mathbf{x},t}^{3/2}$ and flatness $F_{\theta}(y)=\langle\theta^{4}\rangle_{\mathbf{x},t}/\langle\theta^{2}\rangle_{\mathbf{x},t}^{2}$ of temperature, as shown in Fig. 11, where the vertical profiles are averaged over five different realizations of the porous structure at the same porosity. Skewness evaluates the asymmetry of the distribution, whereas flatness quantifies the extent of the distribution’s tails and reveal how extreme values deviate from the mean. Compared to the porous RB cell, the skewness has smaller absolute values near the top (bottom) regions in the canonical RB cell, indicating that the localized cold falling (hot rising) plumes has more profound effects in the canonical RB cell. On the other hand, from Figs. 8(a) and 8(b), we can also observe that the temperature PDF is more symmetric for the porous convection at the heights of $0.25H$ and $0.75H$, which aligns with the lower skewness values seen for temperature in porous convection. In both canonical and porous RB cells, the skewness values are around zero at midheight of $y=0.5H$ [see Figs. 11(a) and 11(c)], indicating almost equal number of hot and cold plumes flow through the midheight. We can also find that the flatness has much smaller values in the porous cell than that in the canonical RB cell, as shown in Figs. 11(b) and 11(d). Specifically, there is a shift towards a Gaussian distribution in the temperature PDF of porous convection as compared to canonical RB convection, because solid porous matrix impedes both hot and cold thermal plumes, there are fewer fluids with temperature that deviate from the bulk temperature. In the canonical RB convection, the differences in skewness and flatness at two Prandtl numbers can be attributed to the shape of large-scale circulation. The LSC is in the form of a tilted ellipse at $Pr=5.3$, occupying a diagonal position within the convection cell, with two secondary corner rolls sited along the opposing diagonal. On the other hand, at $Pr=0.7$, the LSC appears to be circular, with four secondary corner vortices present. This elliptical arrangement of the LSC at $Pr=5.3$ results in asymmetric hot (or cold) plumes falling back to the bottom (or top) plate, which is associated with countergradient heat transfer, thereby creating a more asymmetric temperature PDF, increasing the skewness. At the heights of countergradient heat transfer, there is an abundance of cold and hot plumes with temperatures deviating from the bulk temperature, which are found along the edges of the corner rolls, leading to the peaks in the temperature flatness profile observed at heights of $0.2H$ and $0.8H$.

To evaluate the spatial and temporal distributions of thermal plumes, we adopt the criteria similar to those used in Ref. [53, 54, 25], specifically

Here, $c$ is an empirical constant, for which a value of $c=1$ is chosen. This criterion assumes that plumes occur in regions of local temperature extremes (either maximum or minimum), and in areas where local convective heat flux is larger than the spatial and temporal averaged one. The applicability of this empirical criterion in accurately extracting plume structures within both canonical and porous convection is evident from Fig. 12.

We calculate the time-averaged plume area within the cell, and plot the plume areas as functions of porosity. From Figs. 13 (a) and 13(b), we can see that with the decrease of porosity, plume areas generally increase under both Prandtl number conditions. Additionally, we calculate the plume area along a horizontal line in the fluid phase, and we plot the averaged vertical profile near the bottom wall in Figs. 13(c) and 13(d). These profiles were averaged over five different realizations of the porous structure with the same porosity. Just above the thermal boundary layers ($y\gtrapprox 0.01H$), we observed more hot plumes within the porous convection cell compared to the canonical cell at the same height. This observation serves as another evidence indicating that more hot plumes penetrate the pore throat after detaching from the thermal boundary layers.

A typical snapshot of an instantaneous logarithmic thermal energy dissipation rate field in both the porous and the canonical RB cell is shown in Fig. 14. Overall, we can observe intense thermal energy dissipations occur near the top and bottom boundary layers, where falling cold plumes or rising hot plumes detach from the boundary layers. Besides, we can also observe intense thermal energy dissipation in the bulk region of the porous RB cell, which is absent in the canonical RB cell. The reason behind this observation lies in the permeability of the porous media. In the porous RB cell, plumes can penetrate through the pore throat, leading to the formation of thermal plumes associated with large amplitudes of thermal energy dissipation rates. It should be noted that we assume the solid porous matrix is impermeable to heat flux, thus the thermal energy dissipation rate $\varepsilon_{T}(\mathbf{x})$ is zero at the location of the solid porous matrix, leading to different distributions of $\varepsilon_{T}(\mathbf{x})$ compared to previous study [27].

We further show the time-averaged logarithmic thermal energy dissipation rate field in Fig. 15. We can see that the contribution of thermal plumes to thermal energy dissipation is filtered out in both canonical and porous RB convection cells. In the time-averaged field, we only observe intense thermal energy dissipation occurs near the top and bottom walls, as well as the edge of LSC, where strong temperature gradients exist. Particularly, in the porous cell, we did not observe a significant increase in the thermal energy dissipation rate at the fluid-solid interface of the porous matrix, in contrast to the previous study [27].

We plot the PDFs of thermal energy dissipation rates $\varepsilon_{T}(\mathbf{x}_{f},t)$ obtained over the fluid phase in the cell over time, further normalized by their root-mean-square (rms) values, as shown in Figs. 16(a) and 16(b). Compared with the canonical RB convection, the PDF tails of porous RB convection are more extended, indicating an increasing degree of small-scale intermittency in the thermal energy dissipation field due to the presence of the porous solid matrix. In the figures, we also compared the PDFs of thermal energy dissipation rates with the solid porous matrix being either permeable or impermeable to heat flux. We observed that when the porous matrix is impermeable to heat flux (represented by the red squares), the tails of the PDFs are longer. However, in both scenarios, the tails of the PDFs for thermal energy dissipation rate exceed those found in canonical RB convection. This observation suggests that while the solid porous matrix generally enhance small-scale intermittency, the thermal physical property of the porous matrix also influence the behavior of small-scale intermittency. The PDF of the scalar dissipation rate plays a crucial role in describing turbulent isothermal and reacting flows. It is common to use a log-normal PDF to characterize the distribution of dissipation rate values [55]. We further check whether the thermal energy dissipation fields in the porous RB convection follow a log-normal distribution or a non-log-normal distribution as observed in previous studies of canonical RB convection [25, 24]. In Figs. 16(c) and 16(d), we plot the PDFs of the normalized logarithmic thermal energy dissipation rate $\left(\log_{10}\varepsilon_{T}-\mu_{\log_{10}\varepsilon_{T}}\right)/\sigma_{\log_{10}\varepsilon_{T}}$. We can observe clear departures from log normality of the thermal energy dissipation field for both canonical and porous RB convections, as a result of intermittent local dissipation. Thus, we conjecture that the non-log-normal distribution for thermal energy dissipation rate is universal for buoyancy-driven turbulent convection, even in the presence of complex flow geometry.