# Experimental velocity and temperature measurements for natural convection in a highly porous medium

Nicolò R Sgreva, Christel Métivier, Anthony Teixeira, Tien Dung Le, Sébastien Leclerc

► To cite this version:

Nicolò R Sgreva, Christel Métivier, Anthony Teixeira, Tien Dung Le, Sébastien Leclerc. Experimental velocity and temperature measurements for natural convection in a highly porous medium. International Journal of Thermal Sciences, 2024, 205, pp.109257. 10.1016/j.ijthermalsci.2024.109257. hal-04286807v2

HAL Id: hal-04286807

https://hal.science/hal-04286807v2

Submitted on 18 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Experimental velocity and temperature measurements for natural convection in a highly porous medium

Nicolò R. Sgreva*, Christel Métivier, Anthony Teixeira, Tien Dung Le, Sébastien Leclerc

Université de Lorraine, CNRS, LEMTA, 54000, Nancy, France

## ARTICLE INFO

Keywords:

Solid foam

Natural convection

Magnetic resonance imaging (MRI)

Porous media

## ABSTRACT

We present natural convection experiments at low-to-moderate porous media Rayleigh number $Ra$ in a fluid-saturated structured plastic foam with porosity of 0.83 heated from below. This is motivated by the lack of detailed experiments on thermal convection in these systems, although they are of great interest for many engineering applications. Our experiments were carried out using Magnetic Resonance Imaging (MRI), which allowed us to obtain quantitative measurements of velocity vectors together with temperature fields in the bulk of the fluid-saturated sample. The results obtained are compared with those for classical porous media (e.g. packed beds of spheres) and with theoretical predictions.

For our system, we found the onset of convection at $Ra_c = 41.9$, with a critical wavenumber $\alpha_c \approx \pi$. For $Ra < Ra_c$, velocities always remain below the accuracy of the measurement and temperature fields remain homogeneous and linear in the vertical direction. For $Ra_c \le Ra < 308$, the flow has a stable convective patterns in forms of axisymmetric, straight or deformed rolls. At about $Ra = 308$, the flow enters in a transitional regime in which the convective rolls begin to deform. At the highest $Ra$ tested, the convective pattern takes the form of cross-rolls. The intensity of convection is estimated through the large-scale mean squared velocity $U^2$. For $Ra > 70$, we obtain $U \sim Ra^{1/2}$.

The heat transfer is investigated by evaluating the Nusselt number from temperature fields. We found that the local $Nu$ reflects the flow structure: higher values of local $Nu$ correspond to cold down-welling regions, and smaller $Nu$ to warm up-welling regions. Finally, we found that the averaged Nusselt number for the whole system increases linearly with $Ra$ with a slope equal to 2.

## 1. Introduction

Thermal convection in porous media is a key aspect in many natural phenomena and engineering applications. The former include convective flow of magma in partially crystallized reservoir [1] and groundwater circulation [2]. While in industry, examples include insulation systems [3] and thermal energy storage systems [4]. In the case of thermal energy storage systems, porous media can be combined with phase change materials to form composite systems in which heat is stored and released more efficiently. In this case, the most commonly used macroscopic supports are highly porous open-cell solid foams [5]. In fact, these composite systems can store a large amount of phase change material, have a large surface area to volume ratio, and improve the melting rate [6]. These are made of different materials and typically have a randomly organized solid structure or are characterized by unit cells with a specific geometry [7]. When these composite systems are subjected to a temperature gradient, heat is transferred by conduction and, depending on the system and boundary conditions, by convection.

The competition between the two heat transfer mechanisms depends on the thermal properties of the materials employed, the ability of the fluid to flow between the pores, and the heating position [8,9].

In these cases, a configuration of thermal convection of interest is the porous medium analogue of the Rayleigh-Bénard convection, where a fluid-saturated porous medium is uniformly heated from below and cooled from above. In this case, the fluid flow regime is characterized by the porous medium Rayleigh number ($Ra$), defined as

$$Ra = Ra_f \times Da = \frac{\rho g \beta \Delta T H K}{\eta \kappa_{eff}}, \quad (1)$$

where $Ra_f$ is the Rayleigh number for a non-porous fluid layer and $Da = K/H^2$ is the Darcy number. In Eq. (1), $\rho$ indicates the density of the fluid, $\beta$ its thermal expansion coefficient, $\eta$ its viscosity; $g$ is the acceleration due to gravity, $H$ the thickness of the layer, $\Delta T$ the temperature difference across this thickness, $K$ the permeability and $\kappa_{eff}$ the effective thermal diffusivity of the system. From linear stability

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: nicolo.sgreva@univ-lorraine.fr (N.R. Sgreva).</footnote>
analysis, when Ra remains smaller than a critical value, i.e. Ra < Ra_c, the conductive regime remains stable, the temperature field is uniform and there is no flow [10,11]. The critical value Ra_c = 4π² has been verified several time for different types of porous media and fluids [12]. Large-scale convection in form of cellular motion develops for Ra > Ra_c. For sufficiently small Ra (i.e., Ra_c < Ra < 300), numerical computations based on spectral method show that convection is stable and the convective pattern is characterized by various types of rolls: counter-rotating rolls, superposition of counter-rotating rolls, and polyhedral cells [13]. For larger Ra, the flow is perturbed and becomes less stable due to boundary-layer instabilities. The background cellular structure is preserved until the transition to high-Ra convection from which rolls are completely destabilized [14]. Most experimental work at low-moderate Ra, however, reports irregular polygons from the onset [15], and only a few experiments have succeeded in reproducing the predicted pattern evolution [e.g.16,17]. The latter were performed with classical porous media (packed beds of beads or ordered grids) with relatively low porosity (~0.4–0.5).

The case of natural convection in highly porous systems was studied by Kathare et al. [18] and Davidson et al. [19] who used water-saturated metal and carbon foams. They performed heat transfer measurements and found that the Nusselt number increases with increasing Ra less than what reported for beads. Moreover, for relatively small Ra, the heat transfer enhancement provided by the metal foam occurs mainly by conduction. Shiina and Hishida [20] studied the effects that an anisotropic effective thermal diffusivity has on the onset of convection in a highly porous system formed by thin wires deployed between two horizontal plates. In this case, the critical Rayleigh number was found to depend on both Da and the degree of thermal diffusion anisotropy.

Despite the importance of natural convection, in-depth experiments in porous media, and especially in highly porous systems, remain rare because they are particularly difficult to carry out. In most of them, an idealized granular medium (e.g., packed beds of spheres) or a Hele-Shaw cell is used to model the porous medium [14]. Although the latter allows optical access, it cannot be used to analyze three-dimensional flows and materials with more complex geometries (e.g., with different pore sizes, different pore paths, etc.). On the other hand, the use of a standard porous matrix limits the use of classical optical devices for flow visualization and generally forces the use of invasive measurement techniques. For this reason, most experimental work focuses on measurements of heat transfer within the system, without examining the velocity and temperature fields of the convective fluid.

In recent years, a few experimental techniques have been proposed and tested to resolve fluid flow in the bulk of a fluid-saturated porous material.

Combining heat transfer measurements and optical measurement techniques, Ataei-Dadavi et al. [21] studied natural convection in a porous medium consisting of a bed of packed coarse spheres. Using fluids and porous media with the same refracting index, they have been able to visualize both velocity and temperature fields. They found two regimes of heat transfer: a reduced regime at low Ra, due to the poor penetration of fluid flow through the pores, and a second regime for higher Ra, in which heat transfer is enhanced through the porous medium due to the increased rate of flow penetration across the pores. However, this experimental approach limits the choice of system, which must consist of transparent materials with a refracting index similar to that of the liquid within it.

A different way to visualize flow though an opaque media is by using Magnetic Resonance Imaging (MRI), provided that non-magnetic or electrically conductive materials are used. In this respect, the first attempt was done by Shattuck et al. [22] and Shattuck et al. [17], who used MRI to measure velocity field and investigate convective patterns at various Ra in a packed beds of spheres heated from below and cooled from above. They found that in disordered media, localized convection begins at Ra < Ra_c due to packing defects that lead to locally larger

porosity. In ordered media convection is characterized by convective rolls up to 5Ra_c, but with smaller wavelength than those predicted by theory. Also, for slightly larger Ra (up to 8Ra_c), they obtained a time-dependent state where upward flow regions move randomly in a larger area of downward flow.

Although MRI can provide quantitative information on velocity and temperature, after this first attempt only a few other authors have used this technique to investigate convection in porous media. Weber and Kimmich [23] studied thermal convection in a random-site percolation network using MRI for velocity and temperature mapping. Skuntz et al. [24] and Skuntz et al. [25] analyzed the fluid flow and the melting front of a non-isothermal packed bed of beads containing a phase change material. And Wang et al. [26] carried out experiments of buoyancy-driven convection in an heterogeneous pack of beads. In all these cases, the porous medium was designed by stacking beads together, strongly limiting the maximum possible porosity of the system (the typical porosity for beds of packed spheres is $\phi \sim 0.40$).

In this work we aim to study natural convection at low-to-moderate Ra in an open-cell solid foam with high porosity ($\phi > 0.80$). This is motivated by the lack of detailed experiments on convection in structured solid foam, although it is an extremely important topic for many engineering applications. For this purpose, we carried out experiments using a Kelvin-cell foam as a solid matrix filled with Newtonian liquids. The fluid-saturated foam is then placed in a cylindrical cavity heated from below and cooled from above. We use MRI to measure velocity vectors and temperature fields of the convective fluid within the solid matrix. The materials used during the experiments and their properties are discussed in Section 2. In Section 3 we detail the experimental setup and we describe the protocols used for MRI measurements. Results for the onset of convection, evolution of convection patterns, intensity of convection and heat transport are reported in Section 4. Here we also compare the results obtained for the convection in solid foam with those reported in literature for more classical and idealized porous media. The paper ends with a summary and some concluding remarks in Section 5.

## 2. Materials and their properties

### 2.1. Fluids

Two different Newtonian fluids were used during our experiments: glycerol (85 wt% solution in water) and a paraffin (hexadecane).

The glycerol aqueous solution was prepared by adding demineralized water to pure glycerol, with the latter used as received from the supplier (Fisher Chemical, with a stated purity of 99%). Its thermophysical properties at a temperature of 25 °C are summarized in Table 1. The density $\rho$ was measured at ambient atmospheric pressure with a DMA 5000M Anton-Paar densimeter. The coefficient of thermal expansion $\beta$ was evaluated from the density measurements from $\beta(T_0) = -\rho^{-1}(T_0)(\partial\rho/\partial T)|_{T_0}$, with T being the temperature and $T_0$ the reference temperature. The viscosity $\eta$ was measured with a stress-controlled AR-G2 rheometer (TA instrument) in rotational mode. The measured values (Table 1) are in good agreement with those expected for the same glycerol aqueous solution, see for instance [27].

Hexadecane was supplied by Sigma-Aldrich and labeled as 99% pure. Samples of hexadecane were used as received, i.e. without any further purification or removing impurities. A complete characterization in terms of thermophysical properties of this paraffin is given in [28,29]. The values we report in Table 1 are those given in these last references. The crystallization temperature of hexadecane is $T_s = 17.3$ °C [28]. In all experiments with this fluid, the applied temperatures were always larger than $T_s$ and the material was always in liquid phase.
<figure><img src="image_3.png" /><figcaption>Fig. 1. 3D model used to print the foam sample. The foam disk has a diameter D = 114 mm and a height of H = 26 mm.</figcaption></figure>

**Table 1**

Material properties at T = 25 °C. ρ is the density, η the dynamic viscosity, λ the thermal conductivity and β the coefficient of thermal expansion. The standard uncertainties u_a for the measured properties are u_a(ρ) = 0.01 kg m⁻³, u_a(η) = 5 × 10⁻⁴ Pa s, and u_a(β) = 2.7 × 10⁻⁶ °C⁻¹ with u_a(T) = 0.01 °C.

<table><thead><tr><th>Material</th><th>ρ (kg m<sup>-3</sup>)</th><th>η (Pa s)</th><th>λ (W m<sup>-1</sup> K<sup>-1</sup>)</th><th>β (°C<sup>-1</sup>)</th></tr></thead><tbody><tr><td>85 wt% glycerol</td><td>1215.92</td><td>0.0818</td><td>0.310<sup>a</sup></td><td>5.066 × 10<sup>-4</sup></td></tr><tr><td>Hexadecane</td><td>770<sup>b</sup></td><td>0.003<sup>b</sup></td><td>0.150<sup>b</sup></td><td>8.861 × 10<sup>-4b</sup></td></tr><tr><td>Plastic foam</td><td>1192<sup>a</sup></td><td>–</td><td>0.192<sup>a</sup></td><td>–</td></tr></tbody></table>

The other material properties are from:

a Sgreva et al. [30].

b Sgreva et al. [28].

## 2.2. The porous medium

A structured solid foam with open cells was used as a porous medium. The foam sample was 3D-printed using a light-sensitive photopolymer resin (Monocure 3D Rapid, Monocure Pty Ltd., Australia). The 3D model (Fig. 1) is made by assembling individual Kelvin-type cells to form a cylindrical disk of diameter D = 114 mm and height H = 26 mm. The Kelvin-type unit cell has a diameter (i.e. the distance between two opposing square faces) of 6.5 mm, and cylindrical struts with an approximate diameter of 1.2 mm. The printed foam sample consists then of four vertically stacked cells.

A sample's porosity of φ = 0.83 ± 0.2 was estimated from both magnetic resonance images and via water (volume) absorption. Density and other relevant thermal properties of the solid plastic are reported in [30]. According to it, at a temperature of 25 °C, the density of the solid plastic is 1192 kg m⁻³, the thermal conductivity is 0.192 W m⁻¹ K⁻¹ (Table 1), and the heat capacity is 1.56 kJ kg⁻¹ K⁻¹.

## 2.3. Effective properties of the system

In order to evaluate Ra from Eq. (1), two properties of the volume sample are required: the effective thermal diffusivity κₐ₇₅ of the whole system (i.e. solid matrix and fluid together) and the permeability K of the porous matrix.

Although κₐ₇₅ depends on the thermal diffusivity of the solid matrix and the saturating fluid, expressing it as a weighted arithmetic average of the two contributions is usually not recommended, especially when the thermal conductivity of the different phases is very different from each other [31].

In our case, the effective thermal diffusivity of the system consisting of the 85% aqueous glycerol solution and the solid foam was previously measured experimentally in [30] and the resulting effective thermal diffusivity is κₐ₇₅ = 9.83 × 10⁻⁸ m² s⁻¹ (Table 2). For the case where the saturating fluid is hexadecane a direct measurement

**Table 2**

Effective thermal properties and permeability of the two systems used for the experiments.

<table><thead><tr><th>Property</th><th>Foam + 85wt% glycerol</th><th>Foam + hexadecane</th></tr></thead><tbody><tr><td>Thermal diffusivity (κ<sub>eff</sub>)</td><td>9.83 × 10<sup>-8</sup> m<sup>2</sup> s<sup>-1</sup></td><td>8.97 × 10<sup>-8</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Volumetric heat capacity (ρc<sub>p</sub>)</td><td>2.95 × 10<sup>6</sup> J m<sup>-3</sup> K<sup>-1</sup></td><td>1.75 × 10<sup>6</sup> J m<sup>-3</sup> K<sup>-1</sup></td></tr><tr><td>Permeability (K)</td><td>2.83 × 10<sup>-7</sup> m<sup>2</sup></td><td>2.83 × 10<sup>-7</sup> m<sup>2</sup></td></tr></tbody></table>

is not available. However, as shown by Smith et al. [31] and Sgreva et al. [30], an estimation of κₐ₇₅ can be obtained through the Hashin–Shtrikman model [32]. The latter is based on the volume fraction of the different system components and their thermal conductivities. By using this equivalent model and the weighted average of the volumetric heat capacity of solid and liquid phase, the effective thermal diffusivity for the hexadecane-saturated foam is κₐ₇₅ = 8.97 × 10⁻⁸ m² s⁻¹ (Table 2).

Because of the different viscosities of the saturating fluids and the different effective thermal diffusivities, the two systems used for the experiments result in two different Prandtl numbers, Pr = η/(ρκₐ₇₅). At around T = 25 °C, we obtain a moderate Prandtl number for the hexadecane-saturated foam (Pr = 43.5) and a larger Prandtl number for the glycerol-saturated foam (Pr = 688.1).

For an ideal porous medium made by packed beds of spheres, the permeability K is defined by the Kozeny–Carman equation. In the case of highly porous systems, such as solid foams, the latter model cannot be used and the permeability must be evaluated by experiments or numerical simulations of pressure drop [33]. These experiments have some complexities, such as the difficulty in establishing the flow regime [34] and the proper definition of the morphological parameters used as characteristic length scale [35], that can lead to incorrect estimates of flow characteristics. In our case, the permeability was numerically estimated by solving the Stokes flow on the digital 3D geometry used to print the foam sample. A cubic portion of dimension 26 × 26 × 26 mm with a voxel size of 0.1 mm³ was cropped from the initial 3D image of 114 × 114 × 26 mm (Fig. 1), and was then segmented and meshed by using AVIZO software. COMSOL Multiphysics software based on the finite element method was used to solve the Stokes equations with the Dirichlet boundary condition for the pressure at two opposite faces of the domain in order to obtain the velocity field. The permeability in the flow direction (the direction perpendicular to the faces where the pressure is imposed) can be calculated from the Darcy's law by averaging the velocity field in this direction. Since the sample is isotropic, simulation in only one direction is needed to estimate the permeability value. The value of permeability obtained in this way is K = 2.83 × 10⁻⁷ m² (Table 2), leading to Da = 4.2 × 10⁻⁴. In the same way, we evaluated the tortuosity of the digital porous medium from the effective diffusion coefficient. This led to a value for the tortuosity of the foam of 1.3.
<figure><img src="image_4.png" /><figcaption>Fig. 2. (a) Scheme of half of the experimental setup. We indicate in yellow the cavity containing the sample. (b) Proton density image of a 2 mm horizontal slice (along the xy-plane) and vertical slice (along the xz-plane) of the fluid-saturated solid foam. The horizontal image is taken at half the height of the cavity (i.e. $H/2$, with $H = 0$ on the lower plate of the cavity), while the vertical image crosses the cavity along its central axis.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Diagram (not to scale) of the whole experimental apparatus. $T_1$ and $T_2$ are the two thermocouples that measure the temperature of the upper water flow outlet and inlet, respectively. $T_4$ and $T_5$ do the same for the lower water flow.</figcaption></figure>

# 3. Setup and methods

## 3.1. Experimental setup and protocols

To investigate natural convection in the fluid-saturated solid foam, we used an experimental setup similar to that used in [36,37]. The setup consists of a cylindrical PMMA cavity containing the saturated foam (Fig. 2a). The cavity has the same dimension as the foam disk, that is, an inner diameter of $D = 114$ mm and a height of $H = 26$ mm (aspect ratio $\Gamma = D/H = 4.4$). The solid foam was inserted into the cavity by ensuring good contacts with the internal wall of the PMMA cavity (i.e. no significant gap between foam and cavity wall), and afterward it was carefully filled with the fluid. Air bubbles were removed by manually shaking the device with the help of two small tubes inserted laterally into the cavity through the PMMA ring. The absence of air bubbles was also verified during the experiments directly from magnetic resonance images.

A schematic diagram of the whole experimental apparatus is reported in Fig. 3. Temperatures above ($T_{up}$) and below ($T_{dw}$) the cavity were imposed by water flows from temperature-controlled water baths. Each temperature was measured with two K-type thermocouples (uncertainty of 1%), one placed at the inlet of the water flow and one at the outlet. In this way, $T_{up}$ and $T_{dw}$ indicate the average of the measured inlet and outlet temperatures above and below the cavity, respectively. Good thermal contacts between water flows and the sample were provided by two 3.1-mm-thick sapphire plates with a thermal conductivity of around $35-40 \text{ W m}^{-1}\text{K}^{-1}$. The small thickness and the large thermal conductivity of these plates lead to negligible temperature variations within them. The device was also insulated with a thick layer of insulating foam added to the outside of the side wall.

We applied the following protocol for the experiments. The fluid-saturated foam was initially kept under isothermal conditions ($T_{up} = T_{dw} = T_0$, with $T_0$ being an initial temperature) for a minimum of 5 h. When thermal equilibrium was reached, to obtain convection a temperature difference of $\Delta T = T_{dw} - T_{up} > 0$ was imposed between top and bottom surfaces. $\Delta T$ was afterwards varied by imposing further increments or decrements. The waiting time between each change of $\Delta T$ varied from a minimum of 5 h to more than 12 h. Temperature and velocity fields within the sample were measured at the steady-state. During the entire procedure, the mean temperature $T_m = (T_{dw} + T_{up})/2$ was always kept constant at around $24 \text{ }^{\circ}\text{C}$.

## 3.2. Magnetic Resonance Imaging (MRI)

MRI experiments were performed with a Bruker Avance Biospec 24/40 spectrometer in a temperature-controlled room at $T_{room} = 24 \text{ }^{\circ}\text{C}$. The spectrometer was equipped with a 200-mm wide gradient coil delivering a 200 mT m$^{-1}$ maximum gradient strength and a 160-mm wide quadrature resonator. The horizontal magnetic field and the proton resonance frequency were 2.34 T and 100.3 MHz, respectively.

MRI images were acquired with a slice thickness of 2 mm in a field of view of 120 × 120 mm. Horizontal images were acquired in the xy-plane at mid-height of the cavity (i.e., at $H/2$), while vertical images were acquired in the xz and yz planes along the center axis of the cavity (Fig. 2b). Because of the relatively large extent of the setup relative to the resonator and the gradient coil, the accuracy of MRI measurements is highest in the central part of the images, while the boundaries of the images are characterized by slight distortions of the acquired signal [37]. Moreover, depending on the protocol used and the
<figure><img src="image_6.png" /><figcaption>Fig. 4. MRT calibration for the foam saturated with the 85 wt% glycerol solution at isothermal conditions ($T = T_{up} = T_{dw}$). The figure displays the transverse components of relaxation time $T_2^{NMR}$ as a function of imposed temperature $T$. Data points are fitted by the second order polynomial fit $T_2^{NMR} = 24.161 + 0.835T + 0.043T^2$.</figcaption></figure>

resolution of the images, the acquisition time for each individual image can be relatively long (several minutes). This affects the type of flow regime that can be properly analyzed; i.e., in this case, quasi-steady flow regimes.

MRI was used for three types of measurements: to identify (i) different phases (i.e., the solid foam and the fluid in it), and to obtain (ii) temperature fields and (iii) velocity vectors in the fluid. To do so, the following MRI protocols were used.

(i) Identification and reconstruction of the solid foam geometry were carried out by acquisition of fast proton density maps using Rapid Acquisition with Relaxation Enhancement (RARE) experiments [38]. The identification of the solid structure is possible because of its very fast transverse relaxation time, which leads to absence of signal. Images were acquired with resolutions of 128 × 128 and 256 × 256 pixels (for a minimum pixel length of 0.47 mm). We used an echo time ($t_E$) of 14 ms, a repetition time ($t_R$) of 1000 ms, and a RARE factor of 8.

(ii) Magnetic Resonance Thermometry (MRT) was used to obtain temperature fields of the saturating fluid [37]. Temperature measurements can be achieved when the transverse component of the magnetic relaxation time ($T_2^{NMR}$), i.e. the time required for the magnetization vector to return to its equilibrium state after being subjected to a perturbation, shows a clear dependence with temperature. This depends on the type of fluid used and, in our case, it was observed only for aqueous glycerol solution, but not for hexadecane. For hexadecane, the relaxation time did not show a strong dependence on temperature, thus not allowing MRT to be used. For 85 wt% glycerol, $T_2^{NMR}$ images were obtained by Multi-Slice Multi-Echo (MSME) sequences [37]. We used 12 consecutive echoes, $t_E = 20$ ms and $t_R = 1000$ ms. Images had a resolution of 128 × 128 pixels and were averaged 16 times. The final temperature fields were obtained after a calibration performed under isothermal conditions, in which the correlation between temperature and $T_2^{NMR}$ was established (Fig. 4). More details on the calibration process are given in Supplementary material. The global uncertainty for MRT measurement evaluated during the calibration process is $u(T) = 0.96$ °C, with a standard uncertainty $u_a(T) = 0.95$ °C and an uncertainty on $T_{up}$ and $T_{dw}$ due to thermocouple measurements of 0.10 °C [37].

(iii) Magnetic Resonance Velocimetry (MRV) was used to obtain velocity fields of the fluid within the porous medium. Velocity images were obtained using similar protocols and the same pulse sequence of that of Noel et al. [36]. We used spin-echo imaging sequences together with flow encoding gradient pulses to acquire images with a 128 × 128 pixel resolution. The gradient pulse duration was $\delta = 4$ ms, the gradient

separation time $\Delta_G = 12$ ms, $t_E = 28$ ms, and $t_R = 800$ ms or $t_R = 1000$ ms for experiments with glycerol or hexadecane, respectively. We used a minimum flow encoding gradient pulse intensity of $G_{min} = 0$ mT m$^{-1}$. Because of the different velocities reached during convection in the two fluids, we used a maximum gradient intensity of $G_{max} = 180$ mT m$^{-1}$ for the 85 wt% glycerol solution and $G_{max} = 175$ mT m$^{-1}$ for hexadecane. The velocity field in the direction of the application of the gradient is obtained from the difference between the image acquired with $G_{max}$ and that with $G_{min}$. Following this protocol, the velocity resolution is $5 \times 10^{-5}$ m s$^{-1}$. The global uncertainty for MRV measurements is then estimated from reference images under isothermal conditions in $u(v) = 6.8 \times 10^{-5}$ m s$^{-1}$, with $v$ being the measured velocity.

Differently from MRT where the relaxation time does not depend on the amount of fluid in each voxel, Partial Volume Effect (PVE) arises during MRV if the solid matrix and the saturating fluid are both present in the same voxel. In this case, the voxel signal depends also on the fluid volume fraction in it. In our case, the latter is estimated from normalized RARE images. The normalization is done between the RARE maximum signal, corresponding to a voxel containing only fluid, and the minimum signal corresponding to a voxel of only solid material. The final velocity image, $I_{velocity}$, is then obtained as follows:

$$I_{velocity}^{final} = (I_{velocity}^{raw} - I_{velocity}^{ref}) \times I_{RARE}^{*} \quad (2)$$

where $I_{velocity}^{raw}$ is the velocity image originally acquired, $I_{velocity}^{ref}$ is a reference velocity image obtained without fluid motion ($\Delta T \le 0$), and $I_{RARE}^{*}$ is the normalized RARE image for the fluid volume fraction. Details on MRV calibration are provided in Supplementary material.

# 4. Experimental results

## 4.1. Velocity and temperature measurements and onset of convection

Fig. 5 shows an example of vertical velocity (Fig. 5a) and temperature field (Fig. 5b) within the sample for steady convection at $Ra = 123$ obtained with MRV and MRT, respectively. Here the fluid saturating the solid foam is the glycerol solution and the imposed temperature difference is $\Delta T = 22.3$ °C. Images were acquired 5 h after setting the temperatures. In the vertical scan of Fig. 5b, we also report the contour of three isotherms (20 °C, 25 °C and 30 °C) to better distinguish temperature distribution. Fig. 5 highlights the good agreement between velocity and temperature fields in both planar and vertical scans. In fact, the upward flow (positive vertical velocity) corresponds well to the warmer regions of the fluid, while the downward part of the convective pattern matches the cooler areas.

The convective pattern obtained in Fig. 5 is characterized by convective rolls. This is similar to the results of Shattuck et al. [17] and Howle et al. [16] for ordered porous media but with a lower porosity ($\phi \le 0.5$) than in our case. Moreover, from Fig. 5 it can be seen that the obtained macroscopic flow structure has a larger length scale than the pore scale. This represents one of the conditions needed for the flow to be described by the Darcy's law. The other is a negligible inertia term relative to viscous drag [39]. In our case, this can be estimated by evaluating the local Reynolds number $Re^*$ based on the permeability of the medium, that is

$$Re^* = \frac{|v| K^{0.5} \rho}{\eta}, \quad (3)$$

where $|v|$ is the local velocity magnitude. For an inertia-free flow, $Re^*$ must be small ($Re^* \ll 1$). This translates in terms of Rayleigh number as [14]:

$$\frac{Ra D a^{0.5}}{Pr} \ll 1. \quad (4)$$

The uncertainties for these non-dimensional numbers are discussed in Appendix A.

Fig. 6a shows an example of $Re^*$ calculated from the velocity field obtained from MRV for $Ra = 131.2$. In this case, $Re^*$ remains much
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** $\Delta T = 22.3$ °C, $Ra = 123.0$, $Da = 4.2 \times 10^{-4}$ and $Pr = 688.1$. (a) Vertical velocity field. (b) Temperature field. In the vertical cross-section of the temperature field we highlight the isotherms at 20 °C, 25 °C and 30 °C.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** (a) Local Reynolds number $Re^*$ for convection at $Ra = 131.2$, $Da = 4.2 \times 10^{-4}$ and $Pr = 688.1$. (b) Conditions needed for the flow to be described by the Darcy's law as a function of $Ra$ for all our experiments.</figcaption></figure>

smaller than 1, which indicates that, under these conditions, the flow is inertia-free and can be described by the Darcy's law. Even considering the relatively large pores and the large porosity of our sample foam, this is valid for the majority of the experiments presented here. Only few experiments conducted with hexadecane are near the limit set by Eq. (4), see Fig. 6b. In the following, the latter are used for qualitative purposes only.

Given the good correspondence between temperature field and velocity field at the same $\Delta T$, we estimate the critical difference of temperature ($\Delta T_c$) above which convection occurs directly from the temperature maps in the horizontal $xy$-plane. In fact, under our experimental conditions, temperature variations near the onset are better resolved than velocities, since the latter have very low values.

For the case where the fluid is the aqueous solution of glycerol, $\Delta T_c$ is found by evaluating the variance $\sigma^2$ of the images acquired at different $Ra$. This procedure is similar to that used by Hébert et al. [40] to detect the onset of convection in the classical Rayleigh-Bénard configuration using the shadowgraph method. In our case, to achieve higher accuracy, the variance is computed within a squared region of

length $50 \pm 6$ mm at the center of the temperature image obtained on the $xy$-plane.

The variance of the horizontal temperature field is plotted versus the imposed $\Delta T$ in Fig. 7. Red symbols refer to experiments conducted by increasing $\Delta T$, and blue symbols to those conducted by decreasing $\Delta T$. In both cases, the variance of the temperature field shows an abrupt increase at a specific critical value of $\Delta T$. In the case of Fig. 7, this occurs at $\Delta T_c = 7.6 \pm 0.1$ °C. Below this temperature, $\sigma^2$ is stable around zero, indicating a very homogeneous and constant temperature field within the estimation area. This condition refers to the case where convection has not developed yet. For $\Delta T \ge \Delta T_c$, the variance of the temperature field increases with the imposed temperature difference. The jump of $\sigma^2$ between $\Delta T = 25.7$ °C and $\Delta T = 27.7$ °C (red symbols in Fig. 7) is caused by a change in the convection pattern. We will discuss this aspect in the next section. The value of $\Delta T_c$ and the fluid properties listed in Table 1 lead to a critical Rayleigh number of $Ra_c = 41.9$, with an uncertainty of $u(Ra_c) = 2.1$. The latter is in very good agreement with the theoretical value of $Ra_c = 4\pi^2$ expected for standard porous media.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Variance $\sigma^2$ of the temperature field acquired on the xy-plane as a function of the imposed $\Delta T$. The fluid is the glycerol aqueous solution. Blue symbols refer to decreasing $\Delta T$, while red symbols refer to increasing $\Delta T$.</figcaption></figure>

## 4.2. Convection patterns

The convective pattern in terms of vertical velocity field for different Ra is shown in Fig. 8. Here the Rayleigh number is expressed as reduced Rayleigh-Darcy number $\epsilon$, that is

$$ \epsilon = \frac{Ra - Ra_c}{Ra_c} \qquad (5) $$

In this case, $\epsilon$ was varied from -0.48 to 1.00 (Pr = 688.1). For $\epsilon = -0.48$ (Fig. 8a), no motion is recognizable, the velocity always remains below the accuracy of the measurement, and the image shows only the magnitude of noise. Under similar conditions, Shattuck et al. [17] show that convection can appear along side-walls if porosity increases locally due to structural defects in the porous medium or due to poor thermal insulation. This is not our case, where for $\epsilon < 0$ the onset of convection was never observed either in the bulk of the cavity or along its walls. This is also confirmed by MRT measures. As pointed out previously, at small Ra, MRT provides better visualizations than MRV because less noisy. The temperature fields obtained at $\epsilon < 0$ show an homogeneous and quasi-constant temperature in the regions of the image not affected by signal distortion during scan acquisition. This is shown in Fig. 9a by the normalized temperature profile $T/T_c$ along a cavity's diameter $d$ (Fig. 9b), with $T/T_c$ being the mean temperature along $d$. For $\epsilon = -0.48$, $T/T_c$ remains relatively flat compared to cases with larger $\epsilon$ (Fig. 9a), especially toward the center of the cavity where the effect of signal distortion is minimal, indicating an homogeneous temperature field.

When $\epsilon$ is raised slightly above 0, convection is well established in the bulk of the fluid layer. For $\epsilon = 0.13$, two thick parallel rolls can be recognizable from both velocity (Fig. 8b) and temperature fields (Fig. 9b). Compared with the previous case, now the temperature profile across the convective pattern becomes more wavy, showing a clear increase in temperature variation along $d$ (Fig. 9b). The wavelength $w_L$ associated with the critical condition of motion can be estimated as the distance between the two rolls (e.g. distance between the local maxima of the temperature field). For $\epsilon = 0.13$ (Fig. 9b), we obtain a normalized wavelength $w_L/H$ of 2.04 ± 0.15, leading to a critical wavenumber $\alpha_c = (2\pi H)/w_L = (0.98 \pm 0.08)\pi$, very close to the value $\alpha_c = \pi$ predicted by linear stability theory [12].

In all the experiments we carried out, we found stable rolls up to $\epsilon = 5.43$, corresponding to $Ra = 270.2$. In this range of $\epsilon$ (i.e. $0 < \epsilon < 5.43$) the convective rolls organize themselves mainly as axisymmetric rolls (Fig. 8c), parallel rolls (Fig. 12a), or somewhere in between where one roll is slightly deformed and follows the shape of the container (Fig. 5). These different patterns developed in the same interval of $\epsilon$

are very likely due to variations of thermal boundary conditions at the sidewalls [41].

In the case of axisymmetric rolls, the velocity field of Fig. 8c suggests the presence of small secondary convective structures perpendicular to the sidewall. To better display them, in Fig. 10a we plot only the positive vertical velocity obtained for $\epsilon = 0.64$. Here, between the main axisymmetric roll and the cavity's sidewall, several narrow, elongated, upward-moving structures surrounded by regions of down-welling are visible. Although images resolution does not allow us to evaluate their exact size and geometry, velocities corresponding to these voxels are large enough to be quantified by MRV. This is done in Fig. 10b where the normalized positive vertical velocity is averaged along concentric circles with different radial coordinate $r$ centered in the center of the image. It can be seen that the maximum positive average velocity is obtained at the axis of the main roll, while toward the central region of the image, almost zero positive velocities are measured. On the other hand, in the region near the side wall, the presence of these upward-moving narrow rolls results in an average positive velocity larger than zero and equal to about 15 percent of the maximum average positive velocity. The same phenomena has been observed for classic Rayleigh-Bénard convection. For instance, Croquette et al. [41] observed unstable axisymmetric patterns characterized by the break down of the roll near to the sidewall to form a set of smaller rolls perpendicular to the cylindrical sidewall. In their case, the entire axisymmetric pattern later vanishes into a structure made of fairly regular rolls that always end perpendicular to the sidewalls. However, in our case the axisymmetric pattern remained stable throughout the entire duration of the experiment (up to a maximum of about 60 h after the imposition of $\Delta T$), during which the acquired images showed no changes in the convective pattern.

In Fig. 11 we show the main convective patterns obtained for another experiment where $\epsilon$ was varied from 0.82 to 2.91 (Pr = 688.1). In this case, images were obtained by increasing the previous temperature step by 2 °C. From $\epsilon = 0.82$ to $\epsilon = 2.98$ we obtained the same pattern. In this case, both the velocity field and the temperature field show a warm up-rising roll at the center of the cavity and several warm plumes near the sidewall (Fig. 11a). At the next temperature rise ($\epsilon = 2.64$), the central roll splits into two ends, resulting in polygonal cell pattern (Fig. 11b). Finally, at $\epsilon = 2.91$ the different warm patches re-connect together to form two main rolls almost parallel each others (Fig. 11c). Interestingly, in this case the two rolls do not reach the sidewall, but no smaller structures are visible between them and the sidewall. The pattern change from Fig. 11b to c results in a sudden change in the variance of the temperature field in the xy-plane. This corresponds to the jump of $\sigma^2$ shown in Fig. 7 (red symbols).

The pattern evolution described above was not obtained when $\epsilon$ was decreased from $\epsilon = 2.91$ to $\epsilon = 0.05$. When $\epsilon$ was decreased, the pattern obtained at the highest $\epsilon$ (i.e. the rolls in Fig. 11c) remained stable until $\Delta T$ was reduced to $\Delta T_c$. This hysteresis in pattern evolution may be related to non-Boussinesq effects. These specific experiments were indeed carried out with glycerol solution as the saturating fluid, and the applied temperature difference corresponding to the polygonal cell pattern observed while increasing $\Delta T$ (Fig. 11b) was $\Delta T = 27.7$ °C. A similar hysteresis, in which both rolls and hexagonal convection pattern can be stable at the same Ra, takes place in Rayleigh-Bénard convection when the temperature dependence of fluid properties, such as thermal expansion coefficient and viscosity, is sufficiently large to result in up-down asymmetric conditions in the fluid layer [42,43].

Beside axisymmetric and deformed rolls, from $\epsilon = 2.95$ we observed also stable parallel straight rolls. An example of them is shown in the velocity field of Fig. 12a for hexadecane (Pr = 43.5, $\epsilon = 4.57$). In this case, the wavenumber is $\alpha = (1.14 \pm 0.10)\pi$, slightly larger than $\alpha_c$. This is in agreement with the conditions for stable two-dimensional parallel rolls predicted by theory [44]. In this case, rolls remained stable until $\epsilon \approx 6.32$, corresponding to $Ra \approx 308.0$ and $\Delta T = 1.68$ °C. After that, rolls begin to deform leading to a transitional regime (Fig. 12b). The Ra we
<figure><img src="image_10.png" /><figcaption>Fig. 8. Convective pattern for different $\epsilon$ and $Pr = 688.1$. Here we show the vertical velocity fields for the xy-plane at $H/2$ obtained from MRV. (a) $\epsilon = -0.48$, (b) $\epsilon = 0.13$, (c) $\epsilon = 0.83$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. (a) Temperature profile along the cavity's diagonal $d$ (indicated in (b)) for $\epsilon = -0.48$, $\epsilon = 0.13$, and $\epsilon = 2.91$. The temperature is normalized by the mean temperature $\bar{T}$ along $d$. $d$ is normalized by the cavity diameter $D$. (b) Temperature fields at $\epsilon = 0.13$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Axisymmetric convective pattern. (a) Positive vertical velocity image for half of the xy-plane. $\epsilon = 0.64$, $Pr = 688.1$. The thick black line delimits the wall of the cavity (with $R = D/2$ the cavity's radius). (b) Averaged positive vertical velocity normalized over its maximal value ($v_c^*$) along the radial direction.</figcaption></figure>

obtained for the beginning of this transition matches that predicted by numerical simulations for the transition from a stationary convection regime to a periodic regime characterized by the superposition of counter-rotating rolls. This latter transition, that takes place for $Ra$ in between 300 and 390 [13,45], marks the beginning of thermal boundary-layer instabilities that perturbs (without completely destabilize) the underlying convective structure [14]. When $\epsilon$ is increased further, we finally achieved a clearer cross-roll pattern (Fig. 12c).

### 4.3. Intensity of convection

The intensity of convection for $Ra > Ra_c$ can be evaluated from the velocity image by calculating the squared mean large-scale flow velocity ($U^2$) since $U^2 \sim E_K$, with $E_K$ being the kinetic energy. The large-scale velocity is defined as the surface integral of the magnitude

of the velocity field $v$, that is

$$U^2 = \int_s \|v\|^2 ds, \qquad (6)$$

where $s$ is the area of the velocity image on the $xz$-plane. In our case, $U^2$ is calculated as the sum of the square of the velocity components for each pixel of the velocity image.

Fig. 13 shows the large-scale velocity as a function of $Ra$. Here $U$ is normalized by the characteristic thermal diffusion velocity, $\kappa_{eff}/H$ (with $H$ being the height of the cavity and $\kappa_{eff}$ the effective thermal diffusivity). For $Ra_c < Ra < 70$, a correlation between $U$ and $Ra$ is hard to identify because data are very scattered. This is due to both the difficulty of distinguishing low velocities from background noise and the underestimation of velocity due to the accuracy limit of the MRV. For $Ra > 70$, Fig. 13 shows that the normalized bulk velocity scales
<figure><img src="image_13.png" /><figcaption>Fig. 11. Different patterns for relatively small $\epsilon$. The first row shows the vertical velocity images, while the second row shows the corresponding temperature field (Pr = 688.1). (a) $\epsilon$ = 2.38, (b) $\epsilon$ = 2.64, (c) $\epsilon$ = 2.91.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Convective pattern for hexadecane and different $\epsilon$ (Pr = 43.5). Here we show the vertical velocity fields for the xy-plane at H/2 obtained from MRV. (a) $\epsilon$ = 4.57, (b) $\epsilon$ = 6.32, (c) $\epsilon$ = 25.0.</figcaption></figure>

with Ra as $U \sim Ra^{1/2}$. This is the same scaling obtained by weakly-nonlinear theory for classical Rayleigh-Bénard convection at relatively small Ra [46], and for convection at low Re in porous media made by packed beds of spheres [47].

## 4.4. Heat transport

An aspect of interest in fluid-saturated porous media heated from below is to define the degree to which heat transfer is increased by convection. This can be captured by plotting the Nusselt number (Nu), namely the ratio between the total heat transferred in the system and the heat transferred by conduction alone at the same temperature difference, as a function of $\epsilon$. For $\epsilon < 0$, heat is only transferred by conduction and $Nu = 1$. On the other hand, when convection starts ($\epsilon \ge 0$), the total heat transferred consists of a part due to conduction and a part due to convection leading to $Nu > 1$.

In our experiments with glycerol, Nu can be evaluated from the temperature field obtained by MRT. To do so we used either the vertical cross-section in the xz-plane or the one along yz-plane, depending on the orientation of the convection pattern of each single experiment. In fact, some experiments are characterized by rolls parallel to one of the acquisition planes, limiting their use for evaluating heat transfer. To evaluate Nu we used only those images where at least one convective cell was well-resolved.

Fig. 14 shows an example of the local Nusselt number obtained for $\epsilon$ = 0.62. The local Nu is calculated as

$$Nu = \frac{-\lambda_{\text{eff}} \frac{\partial T}{\partial z}}{-\lambda_{\text{eff}} \frac{\Delta T}{H}}, \qquad (7)$$

where the term $-\lambda_{\text{eff}}(\partial T/\partial z)$ is the total thermal flux, the term $-\lambda_{\text{eff}}(\Delta T/H)$ is the conductive flux through the sample cavity, and $\lambda_{\text{eff}}$ is the effective thermal conductivity of the system. The total thermal flux is obtained from our experiments by evaluating $\partial T/\partial z$ on the first three pixels from the bottom surface of the acquired images. Fig. 14a shows an example of the temperature field on the xz-plane used to calculate Nu. Here the convective rolls form two well-defined ascending warm regions surrounded by descending cold regions. The local Nu (Fig. 14b) calculated for each pixel-column of the temperature image follows this latter trend. Nu takes on relatively higher local values in colder regions, while decreasing as it approaches the warmer regions, with a minimum around the central axes of the latter. The slightly noisy trend for Nu is mainly due to the noise on the temperature field, especially in the first pixels near the bottom boundary. Given the resolution used to acquire these images, the interface between the cavity's bottom surface and the fluid above it is not always well resolved and may lead to some variation in the calculated local Nu.
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Normalized mean large-scale velocity $U$ as a function of Ra. The solid line is the equation $y(x) = ax^{1/2} + b$, with $a = 147.3$ and $b = -774.4$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Pr = 688.1, $\epsilon$ = 0.62 and $\Delta T$ = 13.3 °C. (a) Measured temperature field on the xz-plane during convection. (b) Variation of the local Nusselt number along the normalized cavity's diameter $d/D$.</figcaption></figure>

From Fig. 14b one can evaluate the average Nusselt ($\overline{Nu}$) for a given $\epsilon$ as

$$ \overline{Nu} = \frac{1}{\Sigma} \int_{\Sigma} Nu d\Sigma, \qquad (8) $$

with $\Sigma$ being the cavity's surface. The standard uncertainty is evaluated in $u_a(\overline{Nu}) = 0.1$. For the case in Fig. 14 we obtained $\overline{Nu} = 2.1$. Fig. 15 reports the average Nusselt as a function of $\epsilon$ for all the experiments carried out using glycerol as saturating fluid. As expected, we recovered $\overline{Nu} = 1$ for subcritical conditions and a sharp transition at the onset of bulk convection. For $\epsilon > 0$, $\overline{Nu}$ increases linearly until about $\epsilon = 0.5$ where the slope begins to decrease until it saturates. The difference between the symbols in Fig. 15 reflects the different acquisition planes used to estimate $Nu$. We used images along the xz-plane for experiments developing axisymmetric rolls, and images along the yz-plane for parallel rolls that are perpendicular to the image plane. Results for the axisymmetric case are the most reliable. This is due to the small signal distortion near the sidewall on the xz-plane and the fact that the local $Nu$ can be easily integrated on the cavity's surface $\Sigma$ given the axisymmetric conditions. On the other hand, to avoid the higher lateral distortion, for images of parallel rolls only a portion corresponding to one wavelength of the temperature field on the yz-plane was used. Together with the fact that the convective pattern is not

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Average $\overline{Nu}$ as a function of $\epsilon$ (Pr = 688.1). The different symbols represent the planes on which temperature images were acquired. The solid line shows a slope $S = 2$.</figcaption></figure>

axisymmetric, in this case the local $Nu$ cannot be integrated over the entire surface of the cavity but only along the direction of the profile ($\overline{Nu}_y$). This leads to values of $\overline{Nu}_y$ for parallel rolls that can only reflect an approximation of the real $\overline{Nu}$.

For both parallel and axisymmetric rolls, the results we obtained for the $\overline{Nu}-\epsilon$ correlation in the interval $0 < \epsilon < 0.5$ show a linear trend with a slope $S = d\overline{Nu}/d\epsilon = 2$ (Fig. 15). The same slope is maintained for $\epsilon > 0.5$ in the case of axisymmetric rolls. This is in agreement with the slope $S = 2$ predicted theoretically by Joseph [44] for an unbounded porous layer (aspect ratio $\Gamma \to \infty$) with $Pr Da^{-1} \to \infty$. We recall that the aspect ratio in our case is $\Gamma = 4.4$ and $Pr Da^{-1} = 1.6 \times 10^6$ for glycerol and $Pr Da^{-1} = 1.0 \times 10^5$ for hexadecane.

Because of the accessibility of heat transfer measurements in experiments, this slope has been investigated experimentally in the past for several different systems. For a creeping flow through packed beds of spheres, Elder [48] found a slope of $S = 1$, and Shattuck et al. [17] found $S = 0.73 \pm 0.05$ for aspect ratios $\Gamma = 8-9$. Similar slopes were found by Howle et al. [16], who obtained $0.53 < S < 1.35$ depending on the type of porous medium employed. Again with beads, Close et al. [49] found that $S$ is affected by the aspect ratio of the system, with a decrease in $\Gamma$ leading to an increase in $S$. Beside porous media made by packed spheres, the $\overline{Nu}-\epsilon$ curve for a fluid-saturated highly porous metal foam was investigated by Kathare et al. [18] and Davidson et al. [19]. For $\epsilon > 0$, the value of $S$ they obtained is smaller than what we obtained in our system. This difference may be related to the different geometry and/or different thermal properties of the solid foam. That is, to the different viscous drag between a foam with regular structure and large pores (such as the one used in our experiments) and one with random structure and smaller pores. Or to the departure from the local thermal equilibrium between solid matrix and saturating fluid in the case of a highly conductive metal foam.

## 5. Conclusion

In this work, we carried out experiments of thermal convection at low-to-moderate $Ra$ in a fluid-saturated solid foam heated from below and cooled from above. The solid foam consists of a 3D-printed open-cell structured plastic matrix with a porosity of $\phi = 0.83$ and four unit cells stacked on top of each other. Experimental investigations involving thermal convection in these types of porous media are essential because of their wide use in many heat and mass transfer applications, such as in latent heat thermal energy storage systems.

Natural convection in our system was studied using Magnetic Resonance Imaging (MRI). This technique allowed us to provide for the first
time quantitative measurements of velocity vectors and temperature fields in the convecting fluid. The results obtained were compared with those for classical porous media (e.g. packed beds of spheres with $\phi < 0.5$) and with theoretical predictions.

For our system, we first evaluated the relative importance of the inertia term to the viscous resistance during convection. Even considering the large porous size and high porosity of the foam, for the range of Ra over which the experiments were performed, we obtained $Re^* < 1$ and thus an inertia-free flow that can be effectively described by Darcy's law.

Next, we investigated the critical conditions of motion and the different convection patterns. We found the onset of convection at $Ra_c = 41.9$, with a critical wavenumber of $a_c \approx \pi$, very close to the value predicted by linear stability theory for Rayleigh-Darcy configuration. For $Ra < Ra_c$, we never observed convection. In this case, the velocities always remain below the accuracy of the measurement and temperature fields remain homogeneous and constant. For $Ra_c \le Ra < 308$, the flow organizes in stable rolls, the latter in forms of axisymmetric, straight or deformed rolls. At about $Ra = 308$, the flow enters in a transitional regime in which the convective rolls begin to deform. At the highest $Ra$ tested, the convective pattern takes the form of well-developed cross rolls. The intensity of convection was estimated through the large-scale mean squared velocity $U^2$, since $U^2 \sim E_K$, with $E_K$ as the kinetic energy of the system. For $Ra > 70$, we obtained $U \sim Ra^{1/2}$, as it is for classic Rayleigh-Bénard convection without the porous layer.

Finally, we estimated the Nusselt number from temperature fields. The local $Nu$ reflects the temperature distribution: higher values of local $Nu$ correspond to cold down-welling regions, and smaller $Nu$ to warm up-welling regions. The averaged Nusselt number for the whole system ($\bar{Nu}$) was estimated from the local $Nu$ of images with axisymmetric convective rolls and parallel rolls. We obtained $\bar{Nu} = 1$ for subcritical conditions and a sharp transition at the onset of convection. For experiments with axisymmetrical rolls, when convection begins, $\bar{Nu}$ increases as $\epsilon$ increases with a linear trend with a slope $S = 2$. This is in agreement with the slope predicted by theory, but is different from the slopes obtained in other experiments using different systems. The differences with other experiments reported in the literature may be associated with the different viscous drag that different porous matrices can generate, the departure from local thermal equilibrium if the solid and liquid components have very different thermal properties, and/or the non-homogeneity of the porous medium.

To conclude, in this work, we have successfully coupled MRV and MRT techniques for quantitative and qualitative flow analysis in a highly porous solid system. This approach has allowed us to obtain direct visualization of thermo-convective flows within the bulk of the system, as well as quantitative measurements of velocity and temperature, without perturbing the flow with invasive measurement tools. The results presented here can serve as a base for future studies in opaque porous media, such as the investigation of the impact of different porous geometries on heat and mass transfer and the study of solid-liquid phase change in open-cell foams.

### CRediT authorship contribution statement

**Nicolò R. Sgreva:** Data curation, Formal analysis, Investigation, Validation, Writing – original draft, Writing – review & editing. **Christel Métivier:** Conceptualization, Formal analysis, Funding acquisition, Writing – original draft, Writing – review & editing. **Anthony Teixeira:** Formal analysis, Investigation. **Tien Dung Le:** Investigation, Methodology. **Sébastien Leclerc:** Investigation, Methodology.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgments

This work was founded by the operation STOCK'NRJ, co-financed by the European Regional Development Fund within the framework of the Program FEDER-FSE Lorraine and Massif des Vosges 2019–2023. This work is also part of the project CONVINCES, financially supported by the French National Research Agency (Grant No. ANR-21-CE50-0029).

## Appendix A. Uncertainties for the adimensional numbers

Uncertainties related to material properties measured in this work and those related to MRI techniques are reported in Sections 2.1 and 3.2, respectively. The uncertainty in the non-dimensional numbers is estimated by the error propagation method [50]. For a generic non-dimensional number $R$, this writes

$$u(R) = \sqrt{\sum_{i=1}^{N} \left( \frac{\partial R}{\partial X_i} u(X_i) \right)^2}, \quad (A.1)$$

where $R$ is a function of $N$ independent variables ($X_i$), each with its own uncertainty $u(X_i)$. This leads, for instance, the uncertainty for the Prandtl number (defined as $Pr = \eta/(\rho\kappa_{eff})$) being

$$u(Pr) = \sqrt{\left(\frac{\partial Pr}{\partial \eta} u_a(\eta)\right)^2 + \left(\frac{\partial Pr}{\partial \rho} u_a(\rho)\right)^2 + \left(\frac{\partial Pr}{\partial \kappa_{eff}} u_a(\kappa_{eff})\right)^2}, \quad (A.2)$$

with $u_a(\eta) = 5 \times 10^{-4}$ Pa s, $u_a(\rho) = 0.01$ kg m$^{-3}$, and $u_a(\kappa_{eff}) = 3.5 \times 10^{-9}$ m$^2$ s$^{-1}$ [30]. This results in $u(Pr) = 24.9$ for the glycerol-saturated foam and $u(Pr) = 7.4$ for hexadecane-saturated foam.

Uncertainties for the other non-dimensional numbers are obtained in the same way, i.e. according to Eq. (A.1) and their definitions. For the Rayleigh number, we estimate $u(Ra_c) = 2.1$ at the onset of convective motion for the glycerol-saturated foam. For that we consider, in addition to the uncertainties given above, also $u_a(H) = 0.01$ mm, $u_a(\beta) = 2.6$ °C$^{-1}$ and $u_a(T) = 0.10$ °C for the temperature at the upper and lower cavity's surfaces.

## Appendix B. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.ijthermalsci.2024.109257.

### References

[1] S. Tait, C. Jaupart, Compositional convection in a reactive crystalline mush and melt differentiation, J. Geophys. Res.: Solid Earth 97 (1992) 6735–6756.

[2] H.J. Diersch, O. Kolditz, Variable-density flow and transport in porous media: approaches and challenges, Adv. Water Resour. 25 (2002) 899–944.

[3] D.A.S. Rees, J. Lage, The effect of thermal stratification on natural convection in a vertical porous insulation layer, Int. J. Heat Mass Transfer 40 (1996) 111–121.

[4] F. Agyenim, N. Hewitt, P. Eames, M. Smyth, A review of materials, heat transfer and phase change problem formulation for latent heat thermal energy storage systems (LHTESS), Renew. Sustain. Energy Rev. 14 (2010) 615–628.

[5] S. Zhang, D. Feng, L. Shi, L. Wang, Y. Jin, L. Tian, Z. Li, G. Wang, L. Zhao, Y. Yan, A review of phase change heat transfer in shape-stabilized phase change materials (ss-PCMs) based on porous supports for thermal energy storage, Renew. Sustain. Energy Rev. 135 (2021) 110127.

[6] Y. Yao, H. Wu, Z. Liu, Z. Gao, Pore-scale visualization and measurement of paraffin melting in high porosity open-cell copper foam, Int. J. Therm. Sci. 123 (2018) 73–85.

[7] Z.A. Qureshi, E. Elnajjar, O. Al-Ketan, R.A. Al-Rub, S.B. Al-Omari, Heat transfer performance of a finned metal foam-phase change material (FMF-PCM) system incorporating triply periodic minimal surfaces (TPMS), Int. J. Heat Mass Transfer 170 (2021) 121001.

[8] S. Gasow, A.V. Kuznetsov, Y. Jin, Prediction of pore-scale-property dependent natural convection in porous media at high Rayleigh numbers, Int. J. Therm. Sci. 179 (2022) 107635.
[9] D.J. Krishna, T. Basak, S.K. Das, Natural convection in a non-darcy anisotropic porous cavity with a finite heat source at the bottom wall, Int. J. Therm. Sci. 48 (2009) 1279-1293.

[10] E. Lapwood, Convection of a fluid in a porous medium, in: Mathematical Proceedings of the Cambridge Philosophical Society, Cambridge University Press, 1948, pp. 508-521.

[11] C. Horton, F. Rogers Jr., Convection currents in a porous medium, J. Appl. Phys. 16 (1945) 367-370.

[12] D.A. Nield, A. Bejan, Convection in Porous Media, Vol. 5, Springer, Cham, 2017.

[13] J. Caltagirone, P. Fabrie, Natural convection in a porous medium at high Rayleigh numbers. I: Darcy's model, Eur. J. Mech. B Fluids 8 (1989) 207-227.

[14] D. Hewitt, Vigorous convection in porous media, Proc. R. Soc. A 476 (2020) 20200111.

[15] L. Howle, Convection in ordered and disordered porous layers, in: Transport Phenomena in Porous Media II, Elsevier, 2002, pp. 155-176.

[16] L. Howle, R. Behringer, J. Georgiadis, Convection and flow in porous media. Part 2. Visualization by shadowgraph, J. Fluid Mech. 332 (1997) 247-262.

[17] M. Shattuck, R. Behringer, G. Johnson, J.G. Georgiadis, Convection and flow in porous media. Part 1. Visualization by magnetic resonance imaging, J. Fluid Mech. 332 (1997) 215-245.

[18] V. Kathare, J. Davidson, F. Kulacki, Natural convection in water-saturated metal foam, Int. J. Heat Mass Transfer 51 (2008) 3794-3802.

[19] J.H. Davidson, F. Kulacki, D. Savela, Natural convection in water-saturated reticulated vitreous carbon foam, Int. J. Heat Mass Transfer 52 (2009) 4479-4483.

[20] Y. Shiina, M. Hishida, Critical Rayleigh number of natural convection in high porosity anisotropic horizontal porous layers, Int. J. Heat Mass Transfer 53 (2010) 1507-1513.

[21] I. Ataei-Dadavi, M. Chakkingal, S. Kenjeres, C.R. Kleijn, M.J. Tummers, Flow and heat transfer measurements in natural convection in coarse-grained porous media, Int. J. Heat Mass Transfer 130 (2019) 575-584.

[22] M. Shattuck, R. Behringer, G. Johnson, J. Georgiadis, Onset and stability of convection in porous media: Visualization by magnetic resonance imaging, Phys. Rev. Lett. 75 (1994) 1995.

[23] M. Weber, R. Kimmich, Rayleigh-Bénard percolation transition of thermal convection in porous media: Computational fluid dynamics, NMR velocity mapping, NMR temperature mapping, Phys. Rev. E 66 (2002) 056301.

[24] M.E. Skuntz, D. Perera, J.E. Maneval, J.D. Seymour, R. Anderson, Melt-front propagation and velocity profiles in packed beds of phase-change materials measured by magnetic resonance imaging, Chem. Eng. Sci. 190 (2018) 164-172.

[25] M.E. Skuntz, J.D. Seymour, R. Anderson, Observation of heat transfer due to variable thermophysical properties of sub-, near-and super-critical fluids in porous media by magnetic resonance imaging, Int. Commun. Heat Mass Transfer 128 (2021) 105635.

[26] S. Wang, Z. Cheng, L. Jiang, Y. Song, Y. Liu, Quantitative study of density-driven convection mass transfer in porous media by MRI, J. Hydrol. 594 (2021) 125941.

[27] K. Takamura, H. Fischer, N.R. Morrow, Physical properties of aqueous glycerol solutions, J. Pet. Sci. Eng. 98 (2012) 50-60.

[28] N.R. Sgreva, J. Noel, C. Metivier, P. Marchal, H. Chaynes, M. Isaiev, Y. Jannot, Thermo-physical characterization of hexadecane during the solid/liquid phase change, Thermochim. Acta 710 (2022) 179180.

[29] L. Klochko, J. Noel, N. Sgreva, S. Leclerc, C. Métivier, D. Lacroix, M. Isaiev, Thermophysical properties of n-hexadecane: Combined molecular dynamics and experimental investigations, Int. Commun. Heat Mass Transfer 137 (2022) 106234.

[30] N.R. Sgreva, Y. Jannot, C. Métivier, Effective thermal diffusivity and conductivity of a fluid-saturated solid foam, Rev. Sci. Instrum. 94 (2023).

[31] D.S. Smith, A. Alzina, J. Bourret, B. Nait-Ali, F. Pennec, N. Tessier-Doyen, K. Otsu, H. Matsubara, P. Elser, U.T. Gonzenbach, Thermal conductivity of porous materials, J. Mater. Res. 28 (2013) 2260-2272.

[32] Z. Hashin, S. Shtrikman, A variational approach to the theory of the effective magnetic permeability of multiphase materials, J. Appl. Phys. 33 (1962) 3125-3131.

[33] H. Yang, Y. Li, B. Ma, Y. Zhu, Review and a theoretical approach on pressure drop correlations of flow through open-cell metal foam, Materials 14 (2021) 3153.

[34] N. Dukhan, C.A. Minjeur, A two-permeability approach for assessing flow properties in metal foam, J. Porous Mater. 18 (2011) 417-424.

[35] P. Kumar, F. Topin, State-of-the-art of pressure drop in open-cell porous foams: review of experiments and correlations, J. Fluids Eng. 139 (2017).

[36] J. Noel, C. Metivier, S. Becker, S. Leclerc, Natural convection in phase change material: Experimental study, Int. J. Heat Mass Transfer 183 (2022) 122047.

[37] S. Leclerc, C. Métivier, MRI temperature and velocity measurements in a fluid layer with heat transfer, Exp. Fluids 59 (2018) 1-9.

[38] J. Hennig, A. Nauerth, H. Friedburg, RARE imaging: a fast imaging method for clinical MR, Magn. Reson. Med. 3 (1986) 823-833.

[39] D. Nield, D.D. Joseph, Effects of quadratic drag on convection in a saturated porous medium, Phys. Fluids 28 (1985) 995-997.

[40] F. Hébert, R. Hufschmid, J. Scheel, G. Ahlers, Onset of Rayleigh-Bénard convection in cylindrical containers, Phys. Rev. E 81 (2010) 046318.

[41] V. Croquette, M. Mory, F. Schosseler, Rayleigh-Bénard convective structures in a cylindrical container, J. Physique 44 (1983) 293-301.

[42] F.H. Busse, The stability of finite amplitude cellular convection and its relation to an extremum principle, J. Fluid Mech. 30 (1967) 625-649.

[43] M. Darbouli, C. Métivier, S. Leclerc, C. Nouar, M. Bouteera, D. Stemmelen, Natural convection in shear-thinning fluids: Experimental investigations by MRI, Int. J. Heat Mass Transfer 95 (2016) 742-754.

[44] D.D. Joseph, Stability of Fluid Motions II. Springer, Berlin, Heidelberg, 1976.

[45] S. Kimura, G. Schubert, J.M. Straus, Instabilities of steady, periodic, and quasi-periodic modes of convection in porous media, J. Heat Transfer 109 (1987) 350-355, http://dx.doi.org/10.1115/1.3248087.

[46] A. Schlüter, D. Lortz, F. Busse, On the stability of steady finite amplitude convection, J. Fluid Mech. 23 (1965) 129-144.

[47] E. Palm, J.E. Weber, O. Kvernvold, On steady convection in a porous medium, J. Fluid Mech. 54 (1972) 153-161.

[48] J.W. Elder, Steady free convection in a porous medium heated from below, J. Fluid Mech. 27 (1967) 29-48.

[49] D. Close, J. Symons, R. White, Convective heat transfer in shallow, gas-filled porous media: experimental investigation, Int. J. Heat Mass Transfer 28 (1985) 2371-2378.

[50] R.J. Moffat, Describing the uncertainties in experimental results, Exp. Therm. Fluid Sci. (1988) 13-17.