An explanation to the Nusselt-Rayleigh discrepancy in naturally convected porous media

We model hydrothermal convection using a partial differential equation formed by Darcy velocity and temperature the velocity formulation. Using the Elder problem as a benchmark, we found that the velocity formulation is a valid model of hydrothermal convection. By performing simulations with Rayleigh numbers in the non-oscillatory regime, we show that multiple quasi-steady-state solutions can be one of the reasons that caused the Nusselt-Rayleigh discrepancy found in previous experiments. The results reveal more understandings about the nature of uncertainty of convection modes in porous media.


Introduction
Natural convection in porous media occurs in various geological and industrial settings, such as groundwater, geothermal reservoirs, heat sinks, and thermal energy storage. We use the Rayleigh number (Ra) to characterize the porous media, and the Nusselt number (N u) to measure the quality of convective heat transfer. From a engineering perspective, it is beneficial to correlate the Rayleigh number and the Nusselt number. Cheng (1979) compiled the experimental, analytical, and numerical results of the Nusselt number and the Rayleigh number for convection heat transfer in a porous layer heated from below. The compilation showed widespread Nusselt numbers for a particular Rayleigh number. We refer to this phenomenon as the N u−Ra discrepancy. Experimental and analytical approaches have addressed this problem. Lister (1990) performed experiments using a large porous slab, and showed that lateral thermal dispersion is one of the reasons that caused the N u − Ra discrepancy. Vadasz (2010) and Vadasz and Braester (1992) used analytical techniques to show that boundary and domain imperfections are the causes of the N u − Ra discrepancy. Both authors explain this phenomenon from an experimental point of view that we cannot perform perfect experiments. Therefore, we present numerical simulations, which can be set up as experiments without the imperfections mentioned above.
In a three-dimensional box of saturated porous media, convection can happen in a two or three-dimensional setting (Beck, 1972). Straus and Schubert (1978) found out that two-dimensional flows have a larger Nusselt number compared to three-dimensional flows when Ra ≤ 97 using numerical simulations. The numerical simulations from Holst and Aziz (1972) and Horne (1979) also produced this effect. Therefore, if both two and three-dimensional convection can exist in the same box, then the system can at least have two Nusselt numbers. More promisingly, the simulations of Straus and Schubert (1979) indicate that it is always possible to force either steady two-dimensional or steady three-dimensional convection by proper choice of initial conditions. Govorukhin and Shevchenko (2017) also showed that the selection scenarios of a convection mode strongly depend on the initial temperature distribution of the porous media. The cosymmetry theory developed by Yudovich (1991) showed that there exists an infinite number of stable stationary flow for a fixed Rayleigh number.
Therefore, we hypothesize that the non-uniqueness of convection modes in a box of saturated porous media is one of the reasons that caused the N u − Ra discrepancy. To prove this hypothesis, we perform several numerical simulations using different box sizes and initial conditions of the temperature field. We use the finite element solver -MOOSE Framework (Gaston et al., 2009) to simulate the physics of naturally convected porous media.
Acknowledging that a specific natural or artificial porous media can have multiple heat transfer rates is not too useful for decisions in reservoir engineering. Karani and Huber (2017) used the basin stability analysis (Menck et al., 2013) to quantify the probability of the occurrence of different convection modes in a twodimensional setting. We employ the same method for a three-dimensional setting. This probability of convection modes can provide us more information and obtain the expectation value of the Nusselt number.

Governing Equations
We present the conservation laws that model natural convection in porous media under the following assumptions mentioned by Horne (1979): -The Boussinesq approximation.
-Inertial effects are small, or low Reynolds number.
-Viscosity of the fluid is constant.
Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media -Thermal dispersion is negligible.
-Saturating fluid and porous solid are in thermal equilibrium.
where q = (u, v, w), the non-dimensional Darcy fluxes in x, y and z directions.
T is non-dimensional temperature, the subscript t is the partial derivatives with respect to time, and Ra is the Rayleigh number. The Rayleigh number is defined as where ρ 0 is the reference density, g is the gravity, α is the thermal expansion rate of fluid, K p is the permeability, d is the reference length, ∆T is the temperature difference between the bottom and top boundary, µ is the dynamic viscosity and k m is the overall thermal conductivity. This formulation of velocity and temperature is also used by Florio (2014). The non-dimensional variables are where T 0 is the reference temperature at the bottom boundary, ρ f is the fluid density, p 0 is the reference pressure. The non-dimensional variables are defined with asterisks, and we drop them for convenience throughout this paper. The scalings are chosen such that the square root of the Rayleigh number exists in both momentum and energy conservation equations Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media Since we are using a finite element solver, we require the weak form: The weak form is derived by multiplying the strong form by a test function ω and integrate through the simulation domain. We choose a first order Lagrangian basis function for velocity q and second order Lagrangian basis function for temperature T . The time integration method of the transient solver is the Crank-Nicolson method. We approximate the CFL number on an element as where h min is the minimum length of the element w.r.t x, y and z axes. In general, we pick the maximum CFL number of all elements as the representative CFL number, and we make sure during each timestep, the representative CFL number is around 0.5. When the entropy production is constant with respect to time, we claim that the simulation has reached a steady state (Börsing et al., 2017).

Benchmark Problems and Results
The benchmark results of our code implementation is presented in this subsection. They are essential in this work, but do not serve as our main focus.

The Elder Problem
Elder studied convection caused by localized heating with the Hele-Shaw cell experiment and numerical solutions (Elder, 1967). It is one of the benchmark problems in softwares such as FEFLOW (Trefry and Muffels, 2007), SUTRA (Voss, 1984) and HydroGeoSphere (Brunner and Simmons, 2012) (Simmons and Elder, 2017). Therefore, we benchmark our finite element implementation of the velocity formulation using the Elder problem. The velocity formulation in 2D is The boundary and initial conditions are defined in Figure 1. Using the physical parameters of the Elder problem defined by Graf and Boufadel (2011), we have a Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media non-dimensionalized problem with Ra = 521.3. We use a structured quadrilateral mesh, with the 120 elements in the x direction and 32 elements in the y direction. We choose a uniform timestep ∆t = 0.0001. We compare our solution with the 2, 5, and 10 years simulation results of Graf and Boufadel (2011), which corresponds to 0.01, 0.025 and 0.05 non-dimensionalized time.  Table 1 The physical properties of the porous media. Modified from Graf and Boufadel (2011). The thermal expansion rate is calculated using the density of water as a function of temperature from Thiesen et al. (1900), assuming a linear growth of volume. Beck (1972) derived the preferred cellular mode during the onset of convection of a 3D box full of saturated porous media. We benchmark our code implementations using different boxes with lengths h 1 and depths h 2 . We use the notation of [h 1 , h 2 ] to represent the box dimensions. The boxes we choose are [1.2, 1.2], [2.0, 0.5], and [3.0, 1.0]. The Rayleigh number is set as 42.25, slightly above the critical Rayleigh number, such that the system starts convecting. The initial condition of temperature is the conductive solution with a ±1% perturbation. We use the notation of (m, n) to represent the cellular modes. The results are in Table 2, and they agree with Beck's analytical cellular modes.

Beck's Box
Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media

Quality measures of convective heat transfer
The Nusselt number and entropy production are used to measure the quality of convective heat transfer in our simulations. The Nusselt number is defined as the ratio of convective heat transfer and conductive heat transfer. Due to simplicity and practical reasons, the Nusselt number is widely used in experiments as a quality measure of convection. However, in numerical simulations, we have access to values of temperature and Darcy fluxes in space. Therefore, we can use another quality measurement -entropy production. Entropy production is generally a better assessment of convection, due to its thermodynamic considerations (Herwig, 2016). The Nusselt number is a combination of the quality, given by the entropy production, and the quantity of heat transfer, given by the heat flux involved (Herwig, 2016). Though, in order to revisit the Nusselt-Rayleigh relation and gain new insights from it, our simulations still contain the calculation of the Nusselt Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media where h 1 and h 2 are the length of the box with respect to x and z-direction. This is analogous to the calculations of Hewitt et al. (2014), except that we utilize heat flux on the top boundary. The linear solution of the convective temperature Eq. (12) has zero contribution to the Nusselt number. Therefore, what the Nusselt number actually calculates is the heat flux from the nonlinear part of the temperature T . Bejan (2013) formulates the volumetric rate of the total entropy production in a fluid-saturated porous medium as a result of both irreversible heat transfer (subscript therm) and fluid flow friction (subscript visc) withṠ (J s −1 m −3 K −1 ). The first term on the right-hand side of Eq. (6) represents the entropy production due to heat transfer irreversibility and invokes the rate of heat flow per unit area and unit time, i.e., Fourier's laẇ T top and T bot are the temperature on top and bottom boundaries, respectively. Recall non-dimensional variables Eq. (2) The second term on the right-hand side of Eq. (6) accounts for viscous dissipation effects of the fluidṠ where q and Φ are Darcy flux and the viscous dissipation function, respectively. The second term on the right-hand side of Eq. (9) is only important when the flow tends to behave like non-Darcy flow (Costa, 2006) and can be neglected in our study. In general, the viscous dissipation effects become increasingly important for a heterogeneous medium as preferred fluid flow paths lead to a local increase in the fluid velocity. Take the scaling factor of q in Eq.
(3) and plug in Eq. (9) Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media We integrate entropy production over the computational domain, assuming homogeneous material propertieṡ where V is the volume of the computational domain. Dropping the asterisks, the dimensionless entropy production, or the entropy generation number is We would like to emphasize that the thermal entropy production measures the norm of the linear and nonlinear parts of the temperature gradient, and the Nusselt number measures the flux from nonlinear parts of the temperature from the top boundary. Since the entropy production and the Nusselt number have different meanings in quality measure in convective heat transfer, we calculate both values for all of the simulations.

Problem Formulation
We further extend the entropy production test of a naturally convecting 2D box designed by (Börsing et al., 2017) into a 3D case. The line test and the triangular test is designed in order to provide more insights in the entropy production of naturally convecting 3D boxes. The Rayleigh number is set as 42.25, slightly above the critical Rayleigh number. The characteristic length of the mesh is 0.05. The results of both the transient and steady-state solver is presented and discussed in later chapters. The observations of varying Nusselt number with respect to 2D or 3D cellular modes (Holst and Aziz, 1972) (Straus and Schubert, 1978) (Horne, 1979) can be a reason for the N u − Ra discrepancy. We test this idea by simulating three boxes of different sizes that have 2D or 3D preferred cellular modes during the onset of convection. The long boxes [1.5, 1.0] and [2.3, 0.9] convect with a 2D cellular mode, and the wide box [2.5, 1.5] convect with a 3D one. The three boxes are tested in the regions of 4π 2 < Ra ≤ 196, and the type of convection cells, the Nusselt number, and the entropy production are reported. The characteristic length of the mesh is set to 0.1. The goal is to see how the transition of cellular modes due to increasing Rayleigh number influence the Nusselt number, and whether it can explain the N u − Ra discrepancy.

Simulation Domain and Initial Conditions
We simulate a box of saturated porous media, with a length of h 1 , a depth of h 2 and a height of 1. In general, the initial condition of temperature is a conductive solution adding ±1% perturbation. If we want to force a particular convection mode to occur, we set the initial conditions as Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media  (Beck, 1972). We examine the convection pattern and entropy production of the plotted points using the proposed finite element solver.
Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media where A is the amplitude (Florio, 2014). We set the amplitude such that the entropy production does not exceed 9 256 Ra, which is the analytical upper bound of the N u − Ra relationship (Doering and Constantin, 1998).

Basin Stability Analysis using the Equivalent Entropy Production Initialization
Florio (2014) investigated the probability of convection modes in "critical boxes" using perturbation methods. The critical boxes are boxes with the size that lies on the transition between several convection modes. We investigate this probability using the basin stability analysis similar to Karani and Huber (2017). Instead of exploring the amplitude space from -1 to 1, we consider the symmetry of amplitudes, and we initialize the amplitudes in the space that is bounded by a certain entropy production.
Apply the thermal part of nondimensional entropy production Eq. (11) to T 01 and T 11 , and sum the results Consider the amplitudes A 01 and A 11 as axes in the Cartesian coordinates, the combinations of amplitudes that have the same amount of entropy production NṠ sum is an ellipse where a is the semi-major axis and b is the semi-minor axis. We apply a change of variables such that we can characterize the combination of amplitudes with the same amount of entropy production using the angle θ, and the initial entropy production NṠ sum .
Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media The initial conditions of velocities u, v, w are the same as Eq. (12). The probability of the (0, 1) mode to occur is defined as the area that is characterised by θ = 0 to the separation angle θ = θ * divided by quarter of the elliptical area We divide ξ linearly into 17 parts from 0 to 90 degrees. The entropy productions NṠ sum we choose to initialize are 1.02, 1.04, 1.08, 1.12 and 1.16. Figure 3 contains the nodes of the line test and the triangular test. Figure 5 shows the cellular mode and the entropy production of the line test with Rayleigh number of 42.25. Figure 6 represents the cellular mode and entropy production of the triangular test with Ra = 42.25.

The N u − Ra Relation
We test the Nusselt number of the boxes [1.5, 1.0], [2.3, 0.9] and [2.5, 1.5] with respect to the Rayleigh numbers in the range between 42.25 and 196. Figure 7 and Figure 8 are the results.

The Probability of Mode Occurrence using Basin Stability Analysis
We found out that the separation angle of the box [2 1/4 , 2 1/4 ] between modes (0, 1) and (1, 1) using the basin stability analysis. The separation angle lies between 44.30 • and 49.94 • . We average the two angles and claim the separation angle is 47.12 • . Using Eq. (13), the probability of mode (0, 1) to occur when the initial condition is the superposition of mode (0, 1) and (1, 1) is 46.8 percent. However, consider the symmetry of modes (0, 1) and (1, 0), we can calculate the probability of mode (1, 1) to occur when the initial condition is the superposition of the three modes, which is 36.2 percent. Figure 5 shows that the cellular modes of our simulation match Beck's prediction well. The slight shift of the boundaries of cellular modes change can be attributed to that our simulations are performed on a porous media with Ra = 42.25 > Ra c . We color-coded the cellular modes such that shades of blues and greens are for two-dimensional cellular modes, and gradients of reds and yellows are for threedimensional cellular modes. By only looking at the entropy production of the blues (1,1)

Quality of Convective Heat Transfer in Different Settings of Boxes
(2,1) (3,1) Fig. 5 The line test at Ra=42.25 plotted over Beck's diagram (Beck, 1972). The y-axis of the middle and bottom plots represents the entropy production. The middle plot is solved by the transient solver with the randomly perturbed initial conditions. The bottom plot is solved by the steady-state solver with the initial conditions Eq. (12). The bracketed values (m, n) represent the cellular mode.   Fig. 6 The triangular test at Ra=42.25 plotted over Beck's diagram (Beck, 1972). The bottom plot represents the contour of entropy production using linear interpolation. The bracketed values (m, n) represent the cellular mode.
Focusing on the whole figure including the three-dimensional cellular modes, the observations are: 1. The quality of convective heat transfer of three-dimensional cellular modes are generally worse than those of two-dimensional cellular modes when Ra = 42.25. 2. The entropy production gradually increases from three dimensional cellular modes to two-dimensional cellular modes that have the same total sum of cells m + n.
Non-peer reviewed EarthArXiv preprint submitted to Transport in Porous Media Figure 6 further supports the aforementioned observations. On the boundaries from mode (2, 0) to either (2, 1) or (1, 1), we can also observe a descent of the entropy production. Increasing convection cells do not guarantee a better quality of convective heat transfer. These results give us new insights into the relation between cellular modes and entropy production. We can also view these results as more rigorous benchmarks that compare with theoretical predictions. 4.2 The N u − Ra Relation Figure 8 compared the N u − Ra relation of the three boxes with previous experimental and numerical results. The result in the region of Ra ≤ 100 has some discrepancy, but it does not explain the scattering effect. It can be possible that some experimental settings of porous medium does not allow a certain cellular mode to exist, so the onset of convection happens later when Ra > Ra c = 4π 2 .
Due to the multiple possible steady states of convection cells for Ra > Ra c , the result of 140 < Ra < 196 shows a wide scattering of the Nusselt number from 3.212 to 4.145. Our results agree with Straus and Schubert (1979)'s numerical tests that both two and three-dimensional convection cells could be obtained at 60 ≤ Ra ≤ 150 by perturbing the initial condition of temperature. The results also agree with Borkowska-Pawlak and Kordylewski (1985)'s proof, that continuous transition of pattern flows from two-dimensional to three-dimensional structure is possible (and vice-versa) with Rayleigh number variations.
The hypotheses of box sizes and multiple steady states of convection pattern influencing the N u − Ra relation are proved using the three boxes experiment. Future experiments of the N u − Ra relation should also be aware of the chosen box size and the convection pattern of fluids.

The Probability of Mode Occurrence and its Implications
Using the basin stability analysis with the equivalent entropy production initialization, we can calculate the probability of mode occurrence. In our example of the [2 1/4 , 2 1/4 ] box, the probability of mode (1, 1) to occur is 36.2 percent, which is slightly higher than the modes (0, 1) and (1, 0). From the simulation, we also know that the Nusselt number for mode (0, 1) and (1, 1) is 1.0637 and 1.0546, respectively. Combining with their occurring probability, we can, therefore, calculate the expectation value for the Nusselt number, which is 1.0579. We have provided a method for a new N u − Ra relationship in the non-oscillating region of the Rayleigh number.
The assumption that we can use a straight line to separate the probability space is only for convenience. It is also possible that the line can be of a higher order. As in Figure 9, the straight line does not separate the modes. However, in this example, it already gives a good approximation and demonstrates how we can apply this method to different settings of box sizes and Rayleigh numbers.
We show how the influence of different box sizes and multiple steady states of convection pattern, leads to the discrepancy of N u − Ra relation in the region of low Rayleigh number. We also demonstrate the method of basin stability analysis using equivalent entropy production initialization to study the probability of mode occurrences in a naturally convected porous media. This method can be utilized for further studies of the N u − Ra relationship.

A Methods of Cellular Mode Checking
In the appendix section, we discuss the methods we used for checking cellular modes of a certain convection pattern.

A.1 Counting by Visual Inspection
We can count how many cells a convection pattern has by looking at the temperature profile at y = 0.5. The method is: count how many times of a temperature change and a 0.5 (-) temperature, in either x or z-direction. See Figure 10 for better comprehension of this concept. This process is programmed and used to post process the simulation results in Figure 5, 6, 7 and 8. The problem with this method is that it cannot differentiate between the 3D cellular modes (bottom left of Figure 10) and the linear combination of the 2D cellular modes (bottom right of Figure 10).

A.2 Counting by Fourier Analysis
Our finite element code calculates the temperature in space T (x, y, z) for a given time. The numerical solution consists of the linear solutions Eq. (12), therefore, we can apply the concept of Fourier analysis. We define the convection temperature as the full solution without the conduction solution  (1,1) Fig. 9 The diagram of basin stability analysis. The red line is drawn as y = tan (θ * x), where θ * is the separation angle. The streamlines show how different initial conditions reach the steady-state solution.