Marwan Fahs

Thomas Graf

• Tuong

Vi Tran

Behzad Ataie-Ashtiani

Craig T Simmons

Anis Younes

Study of the Effect of Thermal Dispersion on Internal Natural Convection in Porous Media Using Fourier Series

Keywords: Natural convection, Thermal dispersion, Darcy's law, Fourier series solution, COMSOL multiphysics, Parameter sensitivity analysis, Nusselt number

Natural convection in a porous enclosure in the presence of thermal dispersion is investigated. The Fourier-Galerkin (FG) spectral element method is adapted to solve the coupled equations of Darcy's flow and heat transfer with a full velocity-dependent dispersion tensor, employing the stream function formulation. A sound implementation of the FG method is developed to obtain accurate solutions within affordable computational costs. In the spectral space, the stream function is expressed analytically in terms of temperature, and the spectral system is solved using temperature as the primary unknown. The FG method is compared 1 to finite element solutions obtained using an in-house code (TRACES), OpenGeoSys and COMSOL Multiphysics ® . These comparisons show the high accuracy of the FG solution which avoids numerical artifacts related to time and spatial discretization. Several examples having different dispersion coefficients and Rayleigh numbers are tested to analyze the solution behavior and to gain physical insight into the thermal dispersion processes. The effect of thermal dispersion coefficients on heat transfer and convective flow in a porous square cavity has not been investigated previously. Here, taking advantage of the developed FG solution, a detailed parameter sensitivity analysis is carried out to address this gap. In the presence of thermal dispersion, the Rayleigh number mainly affects the convective velocity and the heat flux to the domain. At high Rayleigh numbers, the temperature distribution is mainly controlled by the longitudinal dispersion coefficient. Longitudinal dispersion flux is important along the adiabatic walls while transverse dispersion dominates the heat flux toward the isothermal walls. Correlations between the average Nusselt number and dispersion 2 coefficients are derived for three Rayleigh number regimes.

Introduction

Natural convection (NC) in saturated porous media is gaining increased attention because of its wide range of applications in petroleum production, geothermal energy, geologic carbon sequestration, nuclear waste disposal, thermal insulation systems, nuclear reactors, oxide fuel cells, solar collectors and biomedical technologies, among many others. Comprehensive reviews on NC in porous media and its applications can be found in [START_REF] Nield | A brief introduction to convection in Porous Media[END_REF], [START_REF] Ataie-Ashtiani | Confusion about "convection[END_REF], [START_REF] Nield | Convection in Porous Media[END_REF]), Vafai (2011, 2015), Vadász (2008) and Ingham and Pop (2005). In the last two decades, numerical modeling has emerged as a key tool to tackle the aforementioned applications. Nowadays, numerical modeling is commonly used to understand NC processes and in several engineering applications such as designing, planning, site investigation, safety report assessment and predictive simulations.

The analysis of many applications of NC in porous media is often based on the constant thermal diffusivity model which neglects thermal dispersion effects. This simplification reduces the computational complexity introduced by the nonlinearity of the thermal dispersion terms.

However, the obvious disadvantage is that the constant thermal diffusivity assumption is not appropriate to describe mixing processes [START_REF] Thiele | Heat dispersion in stationary mixed convection flow about horizontal surfaces in porous media[END_REF][START_REF] Sheremet | Free convection in a porous wavy cavity filled with a nanofluid using Buongiorno's mathematical model with thermal dispersion effect[END_REF][START_REF] Wen | Rayleigh-Darcy convection with hydrodynamic dispersion[END_REF].

Thermal dispersion effects are caused by mechanical mixing of the interstitial fluid at the pore scale. This mixing is mainly related to obstructions introduced by tortuosity of the porous medium. It can also be caused by recirculation generated by low pressure in local regions arising from flow restrictions. A detailed review on thermal dispersion in porous media can be found in [START_REF] Nield | Convection in Porous Media[END_REF]. Thermal dispersion processes are usually modeled as a full velocity-dependent tonsorial quantity, with two components either paral-lel (longitudinal) or orthogonal (transverse) to the flow direction. Several experimental and modeling studies have been conducted to derive the thermal dispersion tensor and/or to estimate the transverse and longitudinal dispersion coefficients depending on the porous media skeleton (e.g., [START_REF] Ozgumus | Effect of pore to throat size ratio on thermal dispersion in porous media[END_REF][START_REF] Degroot | Thermal Dispersion in High-Conductivity Porous Media[END_REF]. [START_REF] Cheng | Thermal dispersion effects in non-Darcian convective flows in a saturated porous medium[END_REF], [START_REF] Hsiao | Natural convection in an inclined porous cavity with variable porosity and thermal dispersion effects[END_REF] and [START_REF] Languri | A combined experimental/numerical approach to study the thermal dispersion in porous media flows[END_REF] have compared numerical simulations to laboratory experiments and have shown that the derived thermal dispersion model provides the best agreement with experimental data. Many studies have highlighted the significant effect of thermal dispersion on natural (e.g., [START_REF] Sheremet | Free convection in a porous wavy cavity filled with a nanofluid using Buongiorno's mathematical model with thermal dispersion effect[END_REF][START_REF] Khaled | Variable porosity and thermal dispersion effects on coupled heat and mass transfer by natural convection from a surface embedded in a non-metallic porous medium[END_REF][START_REF] Abbas | Effect of thermal dispersion on free convection in a fluid saturated porous medium[END_REF][START_REF] Jha | Numerical investigation of transient free convective flow in vertical channel filled with porous material in the presence of thermal dispersion[END_REF], forced (e.g., [START_REF] Kuznetsov | Investigation of the effect of transverse thermal dispersion on forced convection in porous media[END_REF][START_REF] Fahs | Unsteady forced convection heat transfer on a flat plate embedded in the fluid-saturated porous medium with inertia effect and thermal dispersion[END_REF][START_REF] Hooman | Thermal dispersion effects on forced convection in a porous-saturated pipe[END_REF] or mixed convection in non-Darcian porous media (i.e., under a non-Darcian flow regime). Further studies have also shown that strong heat dispersion effects can occur at low pore velocity in the Darcian regime [START_REF] Thiele | Heat dispersion in stationary mixed convection flow about horizontal surfaces in porous media[END_REF][START_REF] Telles | Dispersion in heat and mass transfer natural convection along vertical boundaries in porous media[END_REF][START_REF] Sheremet | Effect of thermal dispersion on transient natural convection in a wavywalled porous cavity filled with a nanofluid: Tiwari and Das' nanofluid model[END_REF];W a n ge ta l .2016; [START_REF] Fajraoui | Analyzing natural convection in porous enclosure with polynomial chaos expansions: effect of thermal dispersion, anisotropic permeability and heterogeneity[END_REF]. Given its significant effects on heat transfer processes and the advancement in computational techniques, the acknowledgement of thermal dispersion effects in the applications related to NC has been recently increased (e.g., [START_REF] Sheremet | Effect of thermal dispersion on transient natural convection in a wavywalled porous cavity filled with a nanofluid: Tiwari and Das' nanofluid model[END_REF][START_REF] Sheremet | Free convection in a porous wavy cavity filled with a nanofluid using Buongiorno's mathematical model with thermal dispersion effect[END_REF][START_REF] Molina-Giraldo | Evaluating the influence of thermal dispersion on temperature plumes from geothermal systems using analytical solutions[END_REF]Özerinçetal.2012). The simulation of thermal dispersion has become available in several software packages such as COMSOL Multiphysics, FLUENT, FEFLOW [START_REF] Diersch | FEFLOW[END_REF] and OpenGeoSys [START_REF] Sachse | OpenGeoSys-Tutorial[END_REF][START_REF] Kolditz | OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media[END_REF], among others.

Thermal dispersion effects have been investigated for problems involving external NC (e.g., [START_REF] Cheng | Thermal dispersion effects in non-Darcian convective flows in a saturated porous medium[END_REF]Plumb1983;[START_REF] Mohammadien | Thermal dispersion-radiation effects on non-darcy natural convection in a fluid saturated porous medium[END_REF][START_REF] El-Hakiem | Thermal dispersion effects on combined convection in non-Newtonian fluids along a nonisothermal vertical plate in a porous medium[END_REF][START_REF] Asbik | The effect of thermal dispersion on free convection film condensation on a vertical plate with a thin porous layer[END_REF]. Thermal dispersion associated to internal NC remains poorly understood, par-3 ticularly relative to studies of external NC or forced convection. Furthermore, despite the increasing interest in the simulation of thermal dispersion processes, there is no reference benchmark solution that can be used to verify the simulation results. The popular problem of NC in square porous cavity is a typical example of internal NC [START_REF] Baïri | A review on natural convection in enclosures for engineering applications. The particular case of the parallelogrammic diode cavity[END_REF]. This problem describes convective flow in a square porous domain subject to differentially heated vertical walls. Due to its simplified boundary conditions and geometry and due to the large number of published numerical solutions based on a variety of numerical methods [START_REF] Fahs | A new benchmark semi-analytical solution for density-driven flow in porous media[END_REF];W a n ge ta l .2018), this problem has been accepted to be a common benchmark for NC codes in the case of constant thermal diffusivity [START_REF] Zhu | Entropy generation due to three-dimensional double-diffusive convection of power-law fluids in heterogeneous porous media[END_REF]. It is also widely used to provide consistent physical understanding of the NC processes in a porous enclosure. To the best of our knowledge, the effect of thermal dispersion on this problem has never been studied. Thus, the main objectives of this work are (1) to develop a high-accuracy solution for the problem of NC in a square porous cavity in the presence of velocity-dependent thermal dispersion, and (2) to use this solution to understand the effect of thermal dispersion on NC processes in a porous enclosure.

The high-accuracy reference solution is derived using the Fourier-Galerkin (FG) spectral element method [START_REF] Peyret | Spectral Methods for Incompressible Viscous Flow[END_REF]. The Galerkin integrals resulting for the dispersion terms are evaluated numerically using a high-order adaptive technique. An efficient approach is developed to solve the spectral system of equations and to deal with the nonlinearity of the convective and thermal dispersion terms. This approach is based on reducing the number of degrees of freedom and on an efficient algorithm for the evaluation of the Fourier series on the integration points. Several test cases, dealing with a variety of parameters, are generated to understand the solution behavior and to provide a set of reference data for code verification. To build more confidence in the FG solution and to highlight its worthiness for code benchmarking, the generated test cases are simulated using an in-house advanced respectively. Vertical and horizontal boundaries are impermeable. Darcy's law is assumed to hold. Flow is considered to be saturated and within the Boussinesq approximation limits.

The porous skeleton is incompressible. Local thermal equilibrium in the porous medium is assumed. Under these assumptions, the governing equations are given by [START_REF] Mahmud | Mixed convection in a square vented enclosure filled with a porous medium[END_REF][START_REF] Rossa | Effects of hydrodynamic dispersion on the stability of buoyancy-driven porous media convection in the presence of first order chemical reaction[END_REF]:

-→ ∇• -→ V 0( 1 ) -→ V - K µ -→ ∇ ( p + ρgz) (2) σ ∂ T ∂t + -→ V • -→ ∇ T - -→ ∇• α m I + α disp -→ ∇ T 0( 3 ) ρ ρ c [1 -β(T -T c )] (4)
where -→ V is velocity vector, p is fluid pressure, T is temperature, K is permeability, µ is fluid viscosity, g is gravitational acceleration, ρ is fluid density, σ is the ratio of the composite material heat capacity to the convective fluid heat capacity, α m is effective thermal diffusivity, I is the identity tensor, α disp is dispersion tensor, ρ c is fluid density at the cold temperature T c and β is thermal expansion coefficient.

Thermal dispersion processes can be described using the Scheidegger model [START_REF] Scheidegger | General theory of dispersion in porous media[END_REF]) that assumes isotropic dispersion processes. In the present work, thermal dispersion is acknowledged via the Bear-Scheidegger model [START_REF] Bear | Introduction to Modeling of Transport Phenomena in Porous Media[END_REF][START_REF] Emami-Meybodi | Dispersion-driven instability of mixed convective flow in porous media[END_REF]. This model suggests that dispersion results in greater mixing in the direction of flow than in the transverse direction. The anisotropic dispersion tensor is defined as:

α disp (α L -α T ) -→ V × -→ V -→ V + α T -→ V I (5)
where α L and α T are the longitudinal and transverse dispersion coefficients, respectively.

The Fourier-Galerkin (FG) Solution

The developed solution is based on the spectral element FG method. Spectral methods combine the exactness of the analytical methods with an important extent of generality in describing the geometry and boundary conditions as with numerical methods (BniLam and Al-Khoury 2017). Spectral methods present a significant interest in accurately solving the equations of flow, mass and heat transfer in porous media [START_REF] Fahs | A new benchmark semi-analytical solution for density-driven flow in porous media[END_REF]Wangetal.2018;[START_REF] Bnilam | A spectral element model for nonhomogeneous heat flow in shallow geothermal systems[END_REF]S h a oe ta l .2015[START_REF] Huddar | Heat transfer in pebble-bed nuclear reactor cores cooled by fluoride salts[END_REF][START_REF] Ameli | Series solutions for saturated-unsaturated flow in multi-layer unconfined aquifers[END_REF][START_REF] Ameli | Semi-analytical 3D solution for assessing radial collector well pumping impacts on groundwater-surface water interaction[END_REF]. Detailed description of the use of spectral methods for solving partial differential equations as well as the implementation of these methods can be found in [START_REF] Kopriva | Implementing Spectral Methods for Partial Differential Equations[END_REF] and [START_REF] Peyret | Spectral Methods for Incompressible Viscous Flow[END_REF].

The solution procedure is described in four steps as outlined below.

Stream Function Formulation and Non-dimensional Equations

The continuity equation can be intrinsically satisfied using the stream function (ϕ), which is defined as:

u ∂ϕ ∂z ; v - ∂ϕ ∂ x ( 6
)
where u and v are the components of the velocity vector ( -→ V ).

The non-dimensional equations can be derived by assuming that buoyancy is the primary force. This can be achieved using the following dimensionless variables:

X x H ; Z z H ; Φ ϕ α m ; θ T -T c T h -T c ; U u • H α m ; V v • H α m (7)
Thus, the stream function formulation of the non-dimensional flow and heat transfer equations (at steady-state conditions) can be written as:

∂ 2 Φ ∂ X 2 + ∂ 2 Φ ∂ Z 2 -Ra ∂θ ∂ X (8) ∂Φ ∂ Z ∂θ ∂ X - ∂Φ ∂ X ∂θ ∂ Z - ∂ 2 θ ∂ X + ∂ 2 θ ∂ Z - ∂ ∂ X A XX disp ∂θ ∂ X + A XZ disp ∂θ ∂ Z - ∂ ∂ Z A XZ disp ∂θ ∂ X + A ZZ disp ∂θ ∂ Z 0( 9 )
where Ra

KHgρ c β(T h -T c ) µα m
is the porous medium Rayleigh number representing the ratio of the buoyancy to diffusion. A XX disp , A XZ disp and A ZZ disp are the coefficients of the non-dimensional dispersion tensor, given as:

A XX disp α xx disp α m A L U 2 √ U 2 + V 2 + R α disp V 2 √ U 2 + V 2 (10) A ZZ disp α zz disp α m A L V 2 √ U 2 + V 2 + R α disp U 2 √ U 2 + V 2 (11) A XZ disp α xz disp α m A L 1 -R α disp UV √ U 2 + V 2 (12)
where A L α L H and R α disp α T α L are non-dimensional longitudinal coefficient and ratio of the transverse dispersion coefficient to the longitudinal dispersion coefficient.

A XX disp , A XZ disp and A ZZ disp can be expressed in terms of the stream function using:

U ∂Φ ∂ Z ; V - ∂Φ ∂ X (13)
Equations ( 8)and( 9) are subject to the following boundary conditions:

X 0 → Φ 0; θ 1 X 1 → Φ 0; θ 0 Z 0, 1 → Φ 0; ∂θ ∂ Z 0( 1 4 )

Change of Variable: Homogeneous Boundary Conditions

Homogeneous boundary conditions are required for the FG method [START_REF] Peyret | Spectral Methods for Incompressible Viscous Flow[END_REF]. This is satisfied for the stream function. For the heat transport equation, as in [START_REF] Fahs | A new benchmark semi-analytical solution for density-driven flow in porous media[END_REF], we applied the following temperature change of variable:

Θ θ + X -1 (15)
Introducing the residuals for flow (R F ) and heat transfer (R H ) equations and using the new variable Θ,Eqs.( 8)and(9) become:

R F ∂ 2 Φ ∂ X 2 + ∂ 2 Φ ∂ Z 2 +Ra ∂Θ ∂ X -Ra 0 (16) R H ∂ 2 Θ ∂ X + ∂ 2 Θ ∂ Z - ∂Φ ∂ Z ∂Θ ∂ X + ∂Φ ∂ X ∂Θ ∂ Z + ∂Φ ∂ Z + A XX disp ∂ 2 Θ ∂ X 2 +2A XZ disp ∂ 2 Θ ∂ X ∂ Z + A ZZ disp ∂ 2 Θ ∂ Z 2 + ∂ A XX disp ∂ X + ∂ A XZ disp ∂ Z ∂Θ ∂ X + ∂ A XZ disp ∂ X + ∂ A ZZ disp ∂ Z ∂Θ ∂ Z - ∂ A XX disp ∂ X - ∂ A XZ disp ∂ Z 0 ( 17
)
The derivatives of the dispersion tensor coefficients with respect to X and Z are calculated as:

∂ A XX disp ∂ X 2 √ U 2 + V 2 A L U ∂U ∂ X + R α disp V ∂ V ∂ X - A XX disp U 2 + V 2 U ∂U ∂ X + V ∂ V ∂ X (18)
∂ A ZZ disp ∂ Z 2 √ U 2 + V 2 A L R α disp U ∂U ∂ Z + V ∂ V ∂ Z - A ZZ disp U 2 + V 2 U ∂U ∂ Z + V ∂ V ∂ Z (19) ∂ A XZ disp ∂ X A L 1 -R α disp √ U 2 + V 2 V ∂U ∂ X + U ∂ V ∂ X - A XZ disp U 2 + V 2 U ∂U ∂ X + V ∂ V ∂ X (20) ∂ A XZ disp ∂ Z A L 1 -R α disp √ U 2 + V 2 V ∂U ∂ Z + U ∂ V ∂ Z - A XZ disp U 2 + V 2 U ∂U ∂ Z + V ∂ V ∂ Z (21)
With the new variable Θ, the thermal boundary conditions become:

X 0, 1 → Θ 0 Z 0, 1 → ∂Θ ∂ Z 0( 2 2 )

The FG Method

The unknowns Φ and Θ now have homogeneous boundary conditions. The problem is also spatially periodic. Thus, Φ and Θ can be approximated by truncated double Fourier series that satisfy the boundary conditions, written as follows:

Φ(X , Z ) N m m 1 N n n 1 B m,n sin(mπ Z )sin(nπ X ) (23) Θ(X , Z ) N r r 0 N s s 1 C r ,s cos(r π Z )sin(sπ X )( 2 4 )
where N m (resp. N r )andN n (resp. N s ) are the number of Fourier series modes (FSMs) in the X and Z directions for the stream function (resp. temperature). B m,n and C r ,s are the Fourier series coefficients for the stream function and temperature, respectively.

The coefficients B m,n and C r ,s can be calculated using a Galerkin-type technique, based on the orthogonality property of the cosine and sine functions. The Fourier series are substituted into Eqs. ( 16)and( 17). Then, the residuals (R F and R T ) are projected into the vector space of the FSMs and set to be zero in the average sense:

R F g,h R F , 4sin(gπ Z )sin(hπ X ) 4 1 0 1 0 R F sin(gπ Z )sin(hπ X )dX dZ 0 g 1 ...Nm; h 1 ...Nn (25) R H g,h R H , 4cos(gπ Z )sin(hπ X ) 4 1 0 1 0 R H cos(gπ Z )sin(hπ X )dX dZ 0, g 1 ...N r ; h 1 ...N s (26)
The Galerkin integrals arising from Eq. ( 25) as well as the ones related to the diffusion and convection terms in Eq. ( 26) [the first five terms in Eq. ( 17)] can be evaluated analytically.

However, because of the nonlinear dependence of the dispersion tensor to velocity components, the integrals related to thermal dispersion cannot be evaluated in a closed form. Thus, we resorted to numerical integration. This yields the following nonlinear spectral system [for more details readers can refer to [START_REF] Peyret | Spectral Methods for Incompressible Viscous Flow[END_REF], [START_REF] Fahs | A new benchmark semi-analytical solution for density-driven flow in porous media[END_REF] and [START_REF] Shao | A high-accurate solution for Darcy-Brinkman double-diffusive convection in saturated porous media[END_REF]]:

R F g,h -π 2 h 2 + g 2 B g,h + Ra π N r r 0 N s s 1 sΛ h,s Λ g,r C r ,s - Ra π 2 Λ g,0 Λ h,0 0 g 1 ...N m , h 1 ...N n (27) R H g,h -ε g π 2 C g,h g 2 + h 2 - π 2 4 N m m 1 N n n 1 N r r 0 N s s 1 s • m • β I g,m,r χ I h,n,s B m,n C r ,s - π 2 4 N m m 1 N n n 1 N r r 0 N s s 1 n.r .β II g,m,r χ II h,n,s B m,n C r ,s + πg Bg,h +4 N p i 1 Wp i F Disp (Xp i , Zp i )cos(gπZp i )sin(hπXp i ) 0, g 0 ...N r , h 1 ...N s (28)
The coefficients of Eqs. ( 27)and( 28) are summarized in "Appendix A". The numerical integrals of the dispersion terms appear as a Riemann sum in the last term of Eq. ( 28). In this sum, N p is the number of integration points, Xp i , Zp i are the coordinates of an integration point (i) and Wp i is the corresponding weight function related to quadrature rule used for numerical integration. F Disp is the function including all dispersion terms, given by:

F Disp (X , Z ) A XX disp ∂ 2 Θ ∂ X 2 +2A XZ disp ∂ 2 Θ ∂ X ∂ Z + A ZZ disp ∂ 2 Θ ∂ Z 2 + ∂ A XX disp ∂ X + ∂ A XZ disp ∂ Z ∂Θ ∂ X + ∂ A XZ disp ∂ X + ∂ A ZZ disp ∂ Z ∂Θ ∂ Z - ∂ A XX disp ∂ X - ∂ A XZ disp ∂ Z (29)
For ultimate accuracy and in order to avoid errors related to the numerical integration, we use an adaptive-order integration technique [START_REF] Mousavi | Efficient adaptive integration of functions with sharp gradients and cusps in n-dimensional parallelepipeds[END_REF]. With this technique, the number of integration points and the order of the Gauss-Kronod rules are automatically updated to reach a prescribed accuracy.

Solving the Spectral System: Numerical Implementation

The FG method provides spectral accuracy for smooth solutions. However, a defect of this method is the Gibbs phenomena [START_REF] Peyret | Spectral Methods for Incompressible Viscous Flow[END_REF][START_REF] Ameli | Series solutions for saturated-unsaturated flow in multi-layer unconfined aquifers[END_REF]) that can lead to instabilities (unphysical oscillations) for sharp solutions, as is the case for NC at high Rayleigh numbers.

In such cases, stable solutions require a large number of FSMs for which the FG method is computationally impractical. [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF] applied the FG method to solve the problem of NC in square porous cavity in the case of constant thermal diffusivity. They developed an efficient implementation of this method that allowed for using a large number of FSMs. Due to this implementation, they obtained solutions for high Rayleigh numbers (up to 10,000).

Here we improve this implementation, and we extend it to deal with a velocity-dependent dispersion tensor. Thus, we reduce the number of degrees of freedom of the nonlinear spectral system by reformulating it with the Fourier coefficients of the shifted temperature (C r ,s )as the primary variables.

A closer look at Eq. ( 27) reveals that coefficients B g,h can be expressed analytically in terms of C r ,s . By substituting B g,h into Eq. ( 28), we can obtain a closed system of equations with C r ,s as the primary variables. Thus, the nonlinear spectral system can be reduced to only the residual vector of the heat transfer equation. To solve this system, we use the efficient nonlinear solver of the IMSL library (http://www.roguewave.com/products-services/imslnumerical-libraries). This solver is based on an advanced Newton's method which is the modified Powell hybrid algorithm. As in [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF], we improve the performance and convergence of the solver by providing the analytical Jacobian matrix and by simplifying the four nested summations using the Kronecker delta properties.

A FORTRAN code is implemented to solve the spectral nonlinear system with the approaches and simplifications developed in the previous section. Some first runs of the code were carried out to verify its performance. These runs revealed that, despite the previous improvements, the code is still very impractical for cases involving velocity-dependent dispersion. For instance, one evaluation of the residual vector would take more than 10 h of CPU time, for N m N r 50 and N n N s 100 and for about 10,000 integration points.

Scrutiny of this matter revealed that the inefficiency of the code is related to the numerical integral procedure used to deal with the dispersion terms. As can be seen in Eq. ( 28), one residual evaluation requires (N r +1) × N s evaluations of the dispersion term (F Disp )onthe integration points. The latter involves the Fourier series for both the stream function and the shifted temperature. This means that, if the standard integration procedure is applied, the Fourier series should be calculated on the integration points (N r +1) × N s times. Thus, the corresponding number of operations is

N p × (N r +1) × N s × [N m × N n + (N r +1) × N s ]
times (i.e., 5 × 10 11 in the previous example). In order to render the code practical, we implement a specific algorithm to reduce the number of evaluations of the Fourier series on the integration points [START_REF] Fahs | The Henry problem: new semianalytical solution for velocity-dependent dispersion[END_REF]). The algorithm proceeds by evaluating F Disp (via the Fourier series and their spatial first-and second-order derivatives) once when evaluating the first component of the heat transport residual vector (R H 0,1 ).The values of F Disp are then stored in a memory array and further used to calculate the other components of the residual vector. With this algorithm, the residual vector can be evaluated with

N p × [N m × N n + (N r +1) × N s ].
For the previous example, this corresponds to 10 8 operations that require only 5 s of CPU time.

Finally, due to symmetrical boundary conditions, we can show analytically that the Fourier series coefficients (B i, j and C i, j ) for which the summation of indices (i + j) is even are equal to zero. We eliminate these coefficients and half the number of degrees of freedom. A common assumption is that λ surf is equal to the equivalent thermal conductivity of the porous medium. The energy balance at the hot wall expresses the continuity of the conductiondispersion flux at X 0 -and the convective flux (calculated via the Newton's cooling law) at X 0 + . The non-dimensional form of the energy balance equation gives:

Nu h λ surf /H 1+A XX disp ∂θ ∂ X X 0 (31)
Knowing that U 0 at the hot wall (impervious boundary conditions), the local Nusselt number can be simplified to:

Nu h λ surf /H 1+A L R α disp |V | ∂θ ∂ X X 0 ( 32
)
This definition of local Nusselt number is similar to the one used in [START_REF] Howle | Natural convection in porous media with anisotropic dispersive thermal conductivity[END_REF] and [START_REF] Sheremet | Effect of thermal dispersion on transient natural convection in a wavywalled porous cavity filled with a nanofluid: Tiwari and Das' nanofluid model[END_REF]. By using the shifted temperature (Θ), Nu becomes:

Nu h λ surf /H 1+A L R α disp |V | ∂Θ ∂ X -1 X 0 ( 33
)
Nu can be evaluated using the Fourier series of Θ and Φ:

Nu h λ surf /H 1+A L R α disp -π N m m 1 N n n 1 nB m,n sin(mπ Z ) π N r r 0 N s s 1 sC r ,s cos(r π Z ) -1 (34)
The average Nusselt number is given by:

Nu 1 0 Nu • dZ (35)
The integral in Eq. ( 35) is evaluated numerically. The temperature distribution and convective fluid flow processes are assessed using the local values of temperature and velocity components at the centers of the top surface (X 0.5; Z 1) and hot wall (X 0; Z 0.5) and at the cavity center within the convective fluid flow. These metrics are calculated using the Fourier series as:

θ top θ (X 0.5; Z 1) N r r 0 N s s 1 (-1) r sin sπ 2 C r ,s +0.5 U top U (X 0.5; Z 1) π N m m 1 N n n 1 mB m,n (-1) m sin(nπ 2); V hot X (X 0; Z 0.5) -π N m m 1 N n n 1 nB m,n sin(mπ 2) (36)

Results and Discussion

Three matters are discussed in this section. We firstly investigate the stability of the FG solution in terms of the FSMs using six test cases. This step is important to confirm that this solution is free of errors. The results of the six test cases are also useful to provide physical insight on the thermal dispersion processes. Secondly, we compare the FG solution against a discontinuous Galerkin finite element (DG) solution obtained using an advanced in-house research code (TRACES) and two standard finite element (FE) solutions based on OGS package and COMSOL Multiphysics ® software. These comparisons are performed to give more confidence in the accuracy of our FG solution and to highlight its worthiness for code benchmarking and validation. Finally, we use the FG solution to develop a complete parameter sensitivity analysis and to understand the role of thermal dispersion coefficients on heat transfer and convective flow processes in a porous square cavity. In our analysis, we consider three levels of Rayleigh numbers (100, 1000 and 5000). As in [START_REF] Wang | Chebyshev spectral element method for natural convection in a porous cavity under local thermal non-equilibrium model[END_REF] and as it is common in the literature, these values are considered to investigate conductiondominated, moderate and convection-dominated cases. We investigate different cases with A L ranging from 0.01 to 1 and R α disp between 0.04 and 1. These ranges of variations are widely used in the literature (e.g., [START_REF] Abarca | Anisotropic dispersive Henry problem[END_REF]), and they are common for field conditions.

Test Cases: Solution Behavior and Understanding Thermal Dispersion Flux

Adequate FSMs should be used to obtain stable FG solutions and to avoid spurious oscillations related to the Gibbs phenomenon. [START_REF] Fahs | A new benchmark semi-analytical solution for density-driven flow in porous media[END_REF][START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF] suggested an efficient technique to select the appropriate FSMs in the case of salinity-driven NC in porous square cavity and seawater intrusion in coastal aquifers, respectively. This technique proceeds by progressively increasing the number of FSMs until reaching a stable solution. The authors showed that the average Sherwood number (solute analogue of the average Nusselt number used in mass transfer) can be used as measure of solution stability. In this work, we adopt the same technique by replacing the Sherwood number by its equivalent in heat transfer, the Nusselt number.

We analyze the behavior of the FG solution via six test-case examples that cover the range of parameters variation. These cases deal with low, intermediate and high Rayleigh numbers (10 2 ,10 3 and 5 × 10 3 , respectively). For each Rayleigh number, we consider two cases involving, respectively, small (A L α L H 0.001) and high (A L 1) longitudinal dispersion coefficients. For all of these cases, and as is common in the literature, we consider R α disp α T α L 0.1 [START_REF] Abarca | Anisotropic dispersive Henry problem[END_REF]).

For a low Rayleigh number (Ra 100), stable solutions for both low and high dispersion cases are obtained using N m N r +1 30 and N n N s 50. As half of the Fourier coefficients are equal to zero (see previous section), the corresponding number of FSMs is 1500. For Ra 1000, in the case of small dispersion coefficient, the buoyancy forces become more important. This sharpens the temperature distribution and concentrates the streamlines in a thin layer toward the walls. The FG solution shows spurious oscillations introduced by the Gibbs phenomenon. The number of FSMs should be increased to remove these oscillations.

In this case, the stable FG solution is obtained using 8000 FSMs (N m N r +1 80 and For all cases, the numerical evaluation of the Galerkin integrals related to the dispersion terms is performed using a specific integration routine of the ISML library based on adaptive higher-order scheme. With this scheme, the domain is subdivided into a certain number of squares and on each square the number of integration points is locally updated to reach the prescribed accuracy. The number of integration points increases with the Rayleigh number.

The solutions for the high dispersion cases involve more integration points than the low dispersion cases. For the low dispersion cases, about 6 × 10 3 points are used for Ra 100, 10 4 points are used for Ra 1000 and 15 × 10 3 points are used for Ra 5000. The corresponding solutions for the high dispersion cases involve about 10 4 ,2× 10 4 and 32 × 10 3 points.

The streamlines and velocity field for the six test cases are plotted in Fig. 2.T h ec o rresponding main isotherms (0.2, 0.4, 0.6 and 0.8) are depicted in Fig. 3. Figures 2 and3 show that for the low dispersion cases, the isotherms and streamlines look much like those obtained using the purely thermal diffusive model [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF]. Closed streamlines and the velocity field show that, due to buoyancy effects, a clockwise-rotating cell is developed in the cavity. A slow-motion circulation zone can be observed at the core of the domain. For the low dispersion case, this zone expands diagonally from the bottom-left corner to the top-right one. The streamlines are highly dense at right-top and left bottom corners, indicating that in these zones the flow velocity reaches its maximum. The isotherms are not vertical but have almost a uniform shape. Due to the circulation of the fluid in place, the hot isotherms tilt upward to the right at the cavity top surface.

For the high dispersion case, the streamlines are less dense compared to the low dispersion case. This indicates that the thermal dispersion process dampens the fluid circulation, which is reasonable because dispersion enhances fluid mixing that dissipates buoyancy forces. This can be particularly observed at the top-left and bottom-right corners where the streamlines are more distant from each other than in the low dispersion case. As a consequence of flow dampened, the zone of slow velocity in the core of the domain expands and becomes less directed toward the diagonal. Despite the fact that mixing processes related to thermal dispersion are dominating, in contrast to the purely diffusive model, the isotherms are not linear. They have an curved shape, which is completely different from those observed with the pure diffusive model. In contrast to the case of low dispersion, the isotherms are not uniformly spaced. To understand the origin of this behavior, we compare the isotherms of the velocity-dependent dispersion model with those obtained using a constant equivalent dispersion tensor. The coefficients of the constant dispersion tensor are calculated as averages of A XX disp , A XZ disp and A ZZ disp over all the domain. For the high dispersion case (Ra 100), the average values are calculated to be 7.8, 0.31 and 7.01, respectively. Figure 4 shows the comparison between the isotherms obtained using both models. It can be clearly seen that with the constant dispersion tensor the isotherms are uniformly spaced as in the pure diffusive model. Thus, this confirms that the zigzag shape of the isotherms is due to the velocity-dependent dispersion processes.

When the Rayleigh number is increased, for the low dispersion case, the streamlines get closer, indicating more intense fluid circulation. The increase in flow velocity shrinks the hydrodynamic boundary layer toward the hot and cold walls. The fluid circulation zone becomes confined toward the walls. The higher Rayleigh number can be seen as an increase in the buoyancy effects at constant diffusion coefficient. Since thermal dispersion is relatively weak, the isotherms become closely controlled by the flow structure. This explains their spiral shape (see Fig. 3). For the high dispersion case, the increase in Rayleigh number increases the flow velocity in the cavity. However, the isotherms are slightly sensitive to the Rayleigh number. This result is in agreement with [START_REF] Wen | Rayleigh-Darcy convection with hydrodynamic dispersion[END_REF]. In that paper, the authors, based on scaling analysis and direct numerical simulations, show that at sufficiently large Rayleigh

Low dispersion case High dispersion case

Ra =100

Ra =1,000 Ra = 5,000 Figure 5a shows the diffusive heat flux and its magnitude as well as the main isotherms. two parts around the central zone (Fig. 5).

In Fig. 6, we distinguish between longitudinal and transverse dispersive fluxes. To understand the behavior of these fluxes, we include the velocity field and the thermal gradient in the same figure. This figure shows that the longitudinal flux that enhances mixing in the direction parallel to the velocity is important in the zones where the velocity is parallel to the thermal gradient (around the top and bottom surfaces, toward the top-right and bottom-left corners). The transverse dispersive flux, which acts in the direction orthogonal to the velocity, is important toward the vertical walls where the velocity is orthogonal to the thermal gradient.

Figure 6a shows that longitudinal dispersion has no impact on heat transfer across the hot wall. This is mainly related to the diffusive flux and to transverse dispersion.

We should mention that, in contrast to standard numerical solution, the results in Figs.

Benchmarking: Comparison with TRACES, OGS and COMSOL

Benchmark reference solutions are valuable to give confidence in the accuracy of the numerical simulations, for the critical assessment of results derived from computer codes and for the improvement in the underlying numerical techniques (e.g., convergence studies). They are valuable to validate newly developed numerical techniques and to compare existing software. Benchmark reference solutions are also helpful for new software users to become familiarized with software interfaces and data processing. For heat transfer applications, reference solutions play a crucial role in testing the simulation procedure as most of the modeling studies are usually performed using a diverse suite of individual computer codes. Thus, benchmark reference solutions are essential to build confidence in the results of these codes and, consequently, on their physical insights and predictive outcomes gained from using them.

In this section, we evaluate the worthiness of the developed FG solutions in benchmarking numerical codes. In fact, an effective benchmark should be unambiguous and representative of physical processes. It is evident that the newly developed FG solution, including velocitydependent dispersion, is more representative of physical processes than existing reference solutions based on the constant thermal diffusivity model. The development of local quantitative metrics provides objective measures of the FG solution and makes it more suitable for benchmarking. However, a common requirement for a useful benchmark is that it allows the examination of the performance of the different numerical techniques implemented in different codes. This section focuses on the evaluation of the ability of the proposed FG solution for benchmarking numerical methods and techniques commonly used in codes. To do so, we compare the FG solution against three numerical codes (TRACES, OGS and COMSOL Multiphysics ® ) based on different numerical techniques. These codes are used. Transient simulations are used in the three codes to avoid convergence issues related to steady-state simulations encountered at high Rayleigh number and in high dispersion cases. Thus, the transient system is allowed to evolve to reach steady-state conditions. TRACES, developed by [START_REF] Hoteit | New two-dimensional slope limiters for discontinuous Galerkin methods on arbitrary meshes[END_REF], is a research code for flow and reactive transport in saturated porous media. This code is based on advanced numerical methods that reduce the numerical artifacts encountered with the standard Eulerian methods while maintaining high performance. It is based on advanced numerical techniques for both spatial discretization and time integration. In TRACES, the flow equations are discretized with the mixed hybrid finite element method. This method is suitable to solve flow in porous media and for coupling mass or heat transfer equation as it provides accurate and intrinsic evaluation of the velocity field. For more details about the mixed hybrid finite element, please refer to the review paper by [START_REF] Younes | Mixed finite elements for solving 2-D diffusion-type equations[END_REF]. The DG method [START_REF] Miller | Numerical simulation of water resources problems: models, methods, and trends[END_REF])isusedforthe discretization of the advection term. This method is suitable for the hyperbolic equation and for solutions involving sharp fronts. It removes oscillations and reduces numerical diffusion [START_REF] Younes | Solving the advection-dispersion equation with discontinuous Galerkin and multipoint flux approximation methods on unstructured meshes[END_REF][START_REF] Nguyen | An implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusion equations[END_REF][START_REF] Dijoux | A projective hybridizable discontinuous Galerkin mixed method for second-order diffusion problems[END_REF]. The diffusion-dispersion term is discretized using the multipoint flux approximation method (MPFA). This method is suitable for anisotropic problems and can be intrinsically coupled with the DG method [START_REF] Younes | Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods[END_REF]. For time integration, TRACES uses high-order methods via an advanced ODE solver based on an adaptive time integration and time stepping scheme. The high-order integration is coupled with the mixed hybrid and DG methods via the method of lines. More details about this code can be found in [START_REF] Younes | Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods[END_REF]. TRACES has been used to simulate NC in porous media in [START_REF] Fahs | A reference benchmark solution for free convection in a square cavity filled with a heterogeneous porous medium[END_REF] and [START_REF] Shao | A new benchmark reference solution for doublediffusive convection in a heterogeneous porous medium[END_REF][START_REF] Fahs | A 3-D semianalytical solution for densitydriven flow in porous media[END_REF]. A structured triangular mesh is used for TRACES simulations. This mesh is obtained by discretizing the domain into several squares and by subdividing each square into four equal triangles.

OGS is a research and scientific framework for the development of numerical methods for the simulation of thermo-hydro-mechanical-chemical processes in porous and fractured media. OGS is a flexible FE framework for solving multi-field problems in porous and fractured media for applications in geoscience and hydrology. In our OGS simulations, we use a non-iterative adaptive time integration scheme based on error estimation [START_REF] Hirthe | Non-iterative adaptive time-stepping scheme with temporal truncation error control for simulating variable-density flow[END_REF]. A uniform square mesh is used for the simulations with OGS.

COMSOL is a simulation software environment for multiphysics simulations. COMSOL is increasingly used in applications related to NC or heat transfer in porous media (i.e., [START_REF] Tan | Closed solutions for transient heat transport in geological media: new development, comparisons, and validations[END_REF][START_REF] Yacine | Soret-driven convection and separation of binary mixtures in a horizontal porous cavity submitted to cross heat fluxes[END_REF]. For the COMSOL simulations, we coupled the "Darcy's law" interface from the "subsurface flow" module to the "heat transfer in porous media" interface from the "heat transfer" module. COMSOL is based on the FE method, but it offers a variety of numerical techniques that can be used to solve the governing equations. In our simulations, we used the default numerical scheme suggested in COMSOL (standard FE is because TRACES is based on the DG method which allows numerical dispersion to be reduced [START_REF] Younes | Solving the advection-dispersion equation with discontinuous Galerkin and multipoint flux approximation methods on unstructured meshes[END_REF][START_REF] Nguyen | An implicit high-order hybridizable discontinuous Galerkin method for linear convection-diffusion equations[END_REF][START_REF] Dijoux | A projective hybridizable discontinuous Galerkin mixed method for second-order diffusion problems[END_REF].

The comparisons developed in this section show the sensitivity of the solutions to the numerical method used for solving the governing equation. The results presented in Table 1 provide high-quality data that is useful for testing numerical codes dealing with NC in porous media. While Fig. 7 provides a qualitative comparison of the solutions, Table 1 allows for a quantitative comparison based on numerical indicators that are more rigorous for benchmarking than visual inspection [START_REF] Prasad | Using quantitative indicators to evaluate results from variable-density groundwater flow models[END_REF].

The buoyancy effects in COMSOL model are implemented by setting up the density as a function of temperature. Thus, the density is dependent on temperature in all terms of the governing equation. This means that the COMSOL model is based on the full nonlinear model without Boussinesq approximation which assumes temperature-dependent density only in the buoyancy term. The good agreement between the COMSOL results and the FG For OGS (mesh of 10 K nodes), inconsistency is related to the adaptive time integration scheme. This error can be removed by imposing a relatively small maximum time step size. For COMSOL (mesh of 2 K nodes), the error appears as unphysical oscillations. These oscillations can be avoided by using a stabilization technique solution provides preliminary insights into the validity of the Boussinesq approximation in the presence of thermal dispersion in NC in a porous enclosure.

Effects of Longitudinal and Transverse Dispersion Coefficients: Correlation of the Nusselt Number with the Dispersion Coefficients

As it can be concluded from the scaled equations, the temperature distribution, heat flux and velocity field are controlled by three non-dimensional parameters: Ra (Rayleigh number),

A L (longitudinal dispersion coefficient) and R α disp (ratio of α T to α L ). Here we use the developed FG solution to examine the role of the dispersion coefficients A L and R α disp on the convective flow and heat transfer processes, for three values of Rayleigh number (100, 1000 and 5000). For each value of Ra, we carried out 48 simulations by uniformly varying A L between 0.01 and 1 and R α disp between 0.04 and 1 (8 values for A L and 6 values for R α disp ). These ranges cover most typical field applications [START_REF] Abarca | Anisotropic dispersive Henry problem[END_REF]). Thus, a total of 144 cases were performed. In all the runs, we used the highest numbers of FSMs (18,000) to be on the safe side when varying the parameters. As a quantity of interest, we consider the temperature (based on square grid 50 × 50 nodes), the Nusselt number and the maximum velocity components (U max and V max ). The 144 simulations developed for the parameter sensitivity analysis have been done efficiently with reduced computational time due to the high computational performance of the FG solution. small (less than 0.15) and a slight decrease for higher A L . This transition regime does not exist in the case of low Rayleigh number and small value of R α disp where U max is continuously decreasing. Figure 9a shows that at constant A L , U max increases with R α disp . Important variations can be observed for the lowest values of R α disp . For the highest values, U max becomes slightly sensitivity to R α disp . V max decreases with A L (Fig. 9b). This confirms that theincreaseinA L enhances the mixing between hot and cold fluids. Such mixing dissipates the buoyancy forces and dampens the fluid circulation. From Fig. 9b, it can be observed that R α disp has a slight effect and in contrast to U max , V max decreases with R α disp . Thus, R α disp has contradictory effects on U max and V max .

Figure 9c shows that, as for the purely thermal diffusive model, the average Nusselt number increases with Ra. This is mainly related to the intensification of the rotating flow associated with the increase of Ra. From Fig. 9c, it can be observed that, at constant A L , the Nusselt number increases with R α disp . This makes sense regarding the expression of .Nu.

(Eq. 32). Despite the fact that the vertical velocity component along the hot wall decreases with increasing A L , this is associated with a decrease of Nu. The increase of A L pushes the 0.9 and 0.1 isotherms toward the vertical wall and as a result increases the temperature gradient which is directly related to the average Nusselt number. The correlation of the average Nusselt number to the Rayleigh number in many configurations of internal NC has been widely investigated in the literature. However, to the best of our knowledge, no correlation study has been done in the presence of thermal dispersion.

Thus, here we present a preliminary correlation study based on the developed FG solution.

In the presence of dispersion, the average Nusselt number is function of Ra, A L and R α disp .

We used the MATLAB curve fitting tool (CFTool) with the results of the 144 runs (presented above) to find this correlation. However, it was quite difficult to find a single correlation for three values of Rayleigh number. Thus, we develop a correlation for each value. The best fittings have been obtained using surface polynomial correlations:

Nu

n A i 1 n R j 1 P i, j A i L R j α ( 37
)
For Ra 100, the results reveal strong correlation (coefficient of determination R 2 1) with a polynomial of third degree (n A + n R 3) (Fig. 12a). When Rayleigh number is increased, the variation of Nu becomes more complex. Higher polynomial degrees should be used in these cases. For Ra 1000 and 5000, good regressions have been obtained using polynomials of degree 4 and 5, respectively (Fig. 12b,c). The polynomial coefficients for

Conclusion

This study presents a comprehensive analysis of the effect of thermal dispersion on internal NC in a square porous enclosure. A vertically oriented porous square cavity with heated vertical walls is investigated. Thermal dispersion is modeled using a velocity-dependent anisotropic dispersion tensor, introducing longitudinal and transverse dispersion coefficients.

A high-accuracy solution is developed using the FG method applied to the stream function formulation of the governing equation. A sound implementation is developed to obtain accurate solutions with optimized CPU time. In the spectral space, the velocity field is expressed analytically in terms of temperature, and the spectral system is solved with the Fourier coefficients of the temperature as the primary unknown. Comparisons against TRACES (in-house

Fig. 1

 1 Fig. 1 Conceptual model of the studied problem including heat boundary conditions. All boundaries are impermeable to fluid flow

  :2019/10/30 Pages: 33 Layout: Small u n c o r r e c t e d p r o o f Study of the Effect of Thermal Dispersion on Internal…

  o r r e c t e d p r o o f

  :2019/10/30 Pages: 33 Layout: Small u n c o r r e c t e d p r o o f Study of the Effect of Thermal Dispersion on Internal…

3. 5

 5 Quantitative Indicators for Code Verification and Benchmarking Prasad and Simmons (2005) highlighted the importance of quantitative indicators for code benchmarking compared to direct comparison of state variables maps due to the complexity of the instability processes and the qualitative nature of visual comparisons. Thus, we propose several objective metrics to assess heat transfer and flow processes which could be useful as quantitative indicators for code validation and benchmarking. The heat transfer processes are assessed using the average Nusselt number. The local Nusselt number (Nu) provides a measure of the convection heat flux occurring at the hot wall. It is defined as follows: Nu h λ surf /H (30) where h is the local convection coefficient and λ surf is the thermal conductivity of the material forming the hot wall surface. 123 Journal: 11242 Article No.: 1356 TYPESET DISK LE CP Disp.:2019/10/30 Pages: 33 Layout: Small u n c o r r e c t e d p r o o f Study of the Effect of Thermal Dispersion on Internal…

  the high dispersion case, the solution is smoothed by the mixing processes induced by thermal dispersion. The Gibbs phenomenon completely disappears and a stable solution can be obtained with the same number of FSMs as in the low Rayleigh number case (N m N r +1 30 and N n N s 50). For Ra 5000, 18,000 FSMs (N m N r +1 120 and N n N s 150) should be used to obtain a stable solution for the low dispersion case. For the high dispersion case, the FG solution is obtained using 3070 FSMs (N m N r 30, N n 60 and N s 140). 123 Journal: 11242 Article No.: 1356 TYPESET DISK LE CP Disp.:2019/10/30 Pages: 33 Layout: Small u n c o r r e c t e d p r o o f Study of the Effect of Thermal Dispersion on Internal…

Fig. 2 Fig. 3

 23 Fig. 2 The streamlines and velocity field for six test examples dealing with Ra 100, 1000 and 5000 in the case of low dispersion (A L 0.001) and high dispersion (A L 1) coefficients (R α disp 0.1 for all cases)

Fig. 4

 4 Fig. 4 Main isotherms (0.2, 0.4, 0.6 and 0.8) for Ra 100, A L 1andR α disp 0.1 (high dispersion case): comparison between the velocity-dependent dispersion model (color flood map) and the constant dispersion model (solid lines)

Figure

  Figure5bdepicts the same patterns for the dispersive heat flux. Figure5shows that, for the case considered, the magnitude of the dispersive flux is one order of magnitude greater than that for the diffusive flux. The mixing between hot and cold fluid is dominated by thermal dispersion. Different behaviors can be observed for both fluxes. The diffusive flux is important around the highest and lowest isotherms in particular toward the bottom-left and top-right corners. Dispersive flux is important along the adiabatic walls. It exhibits its maximum amplitude toward the bottom-left and top-right corners (highest and lowest isotherms). The velocity is relatively weak at the center of the convective cell. Thus, the dispersive flux is almost zero in this zone. The dispersive heat flux coming from the hot wall is separated into

5 and 6 Fig. 5 Fig. 6

 5656 Fig. 5 Spatial distribution of the diffusive heat flux (a) and dispersive heat flux (b)forRa 100, A L 1and R α disp 0.1 (high dispersion case): flux magnitude (color flood map), flux (arrows) and isotherms (contours)

123

  Journal: 11242 Article No.: 1356 TYPESET DISK LE CP Disp.:2019/10/30 Pages: 33 Layout: Small u n c o r r e c t e d p r o o f Study of the Effect of Thermal Dispersion on Internal…

Fig. 8

 8 Fig. 8 Temperature distribution: examples of inconsistent results obtained with OGS (left, Ra 1000) and COMSOL (right, Ra 5000).For OGS (mesh of 10 K nodes), inconsistency is related to the adaptive time integration scheme. This error can be removed by imposing a relatively small maximum time step size. For COMSOL (mesh of 2 K nodes), the error appears as unphysical oscillations. These oscillations can be avoided by using a stabilization technique

Figure 9 Fig. 9

 99 Figure 9 depicts the variations of U max , V max and Nu for the three Rayleigh number regimes. The first observation is that U max and V max both increase with Ra. Thus, as with the purely diffusive model, the increase in Rayleigh number intensifies the rotating convective flow no matter what the range of the dispersion coefficients. This is consistent with the analysis in the previous section based on streamlines. U max is almost insensitive to A L , especially at high Rayleigh numbers. For the high and intermediate Rayleigh numbers, at constant R α disp , the effect of A L on U max indicates two regimes. We have increasing variation when A L is 123 Journal: 11242 Article No.: 1356 TYPESET DISK LE CP Disp.:2019/10/30 Pages: 33 Layout: Small

  Fig. 11 Variations of the main isotherms (0.1, 0.5 and 0.9) for Ra 100: a) with respect to A L at constant R α disp ( 0.4) and b) with respect to R α disp at constant A L ( 0.3). The arrows indicate increasing variation of A L and R α disp

  Fig. 12 Variation of the average Nusselt number with respect to A L and R α disp with the corresponding regression surface for a Ra 100, b 1000 and c 5000

  code), OGS and COMSOL highlight the worthiness of the developed FG solution as a benchmark to validate and compare numerical codes. For low Rayleigh numbers, the FG solutions are in excellent agreement with TRACES, COMSOL and OGS. For high Rayleigh numbers, the FG solution shows high accuracy as it avoids numerical artifacts related to spatial and time discretization, encountered with standard methods. High-quality data based on quantitative indicators, suitable for code benchmarking, are provided. The FG solution is used to obtain 123 Journal: 11242 Article No.: 1356 TYPESET DISK LE CP Disp.:2019/10/30 Pages: 33 Layout: Small

Table 1

 1 The non-dimensional metrics obtained using the FG solution, TRACES (in-house code), OpenGeoSys (OGS) and COMSOL for the low dispersion case (referred to as "L-disp") and the high dispersion case (referred to as "H-disp"). "Er" is the relative error calculated by assuming the FG values' as error-free references

		Code	Nu	Er Nu	θ top	Er θ top	U top	Er U top	V top	Er V top
	Ra 100
	L-Disp	FG	3.14	-	0.86	-	17.67		35.94
		TRACES	3.14	0.0%	0.86	0.0%	17.74	0.39%	35.20	2.05%
		OGS	3.11	0.95%	0.86	0.0%	17.90	1.30%	36.64	1.94%
		COMSOL	3.11	0.95%	0.86	0.0%	17.92	1.41%	36.55	1.69%
	H-Disp	FG	7.98	-	0.59	-	24.40	-	37.11	-
		TRACES	8.01	0.37%	0.59	0.0%	24.32	0.32%	36.57	1.45%
		OGS	8.05	0.87%	0.59	0.0%	24.55	0.61%	37.61	1.91%
		COMSOL	8.01	0.37%	0.59	0.0%	24.61	0.86%	37.82	1.34%
	Ra 1000
	L-Disp	FG	15.51	-	0.94	-	83.26	-	408.39	-
		TRACES	15.72	1.35%	0.94	0.0%	83.43	0.20%	400.33	1.97%
		OGS	15.09	2.70%	0.94	0.0%	84.66	1.68%	417.38	2.20%
		COMSOL	15.12	2.51%	0.94	0.0%	84.69	1.71%	420.92	3.06%
	H-Disp	FG	69.01	-	0.6	-	232.06	-	381.02	-
		TRACES	68.55	0.66%	0.6	0.0%	234.39	1.00%	380.45	0.14%
		OGS	69.14	0.36%	0.6	0.0%	234.04	0.85%	387.01	1.57%
		COMSOL	69.26	0.18%	0.6	0.0%	235.43	1.45%	387.43	1.68%
	Ra 5000
	L-Disp	FG	66.05	-	0.97	-	301.41	-	2173.90	-
		TRACES	62.72	5.04%	0.96	1.03%	297.84	1.18%	2107.88	3.03%
		OGS	58.11	12.00%	0.96	1.03%	305.12	1.23%	2254.12	3.69%
		COMSOL	57.18	13.42%	0.96	1.03%	304.90	1.15%	2236.25	2.86%
	H-Disp	FG	341.95	-	0.61	-	1154.32	-	1910.64	-
		TRACES	337.04	1.43%	0.61	0.0%	1167.61	1.15%	1899.83	0.56%
		OGS	339.11	0.83%	0.61	0.0%	1159.91	1.28%	1914.04	0.17%
		COMSOL	343.04	0.31%	0.61	0.0%	1169.21	0.48%	1924.10	0.70%

Author Proof u n c o r r e c t e d p r o o f

Study of the Effect of Thermal Dispersion on Internal… method for spatial discretization, backward differentiation formula for time integration and fully coupled approach). An unstructured triangular mesh is used in COMSOL simulations.

The same test cases as in the previous section are simulated using TRACES, COMSOL and OGS. The physical parameters used in the three codes are summarized in "Appendix B".

Computational meshes involving an equivalent number of nodes are used in the three codes (about 10 K nodes).

Figure 7 illustrates the isotherms obtained using the FG method and the three codes (TRACES, COMSOL and OGS). The corresponding non-dimensional metrics (Nu, θ top , U top and V top ) are summarized in Table 1. Also included in this table are the numerical errors for 4 the three codes when considering the FG solution as the reference (error-free) solution. For Ra 100 and in both low and high dispersion cases, Fig. 7 shows excellent agreement between all solutions regarding the isotherms. With errors less than 2%, Table 1 confirms the excellent agreement between the different solutions regarding the non-dimensional metrics. This result gives confidence in the accuracy of (1) the code developed for the FG solution, (2) the procedure used to select the FSMs, and (3) the technique used to include velocity-dependent dispersion (i.e., the numerical evaluation of Galerkin integrals related to the dispersion terms).

The cases with Ra 100 are simple from a computational point of view. The three codes achieve acceptable solutions even with relatively coarse meshes.

For higher Rayleigh (Ra 1000 or 5000), with a relatively coarse mesh, OGS exhibits spurious/unphysical oscillations in the low dispersion cases. The same problem can be observed with COMSOL (Fig. 8). For equivalent computational meshes (same number of nodes), the oscillations with the square grid (in OGS) are larger than those with the triangular grid (in COMSOL). These oscillations are related to the standard FE method. They do not appear in TRACES as this is based on the DG method. To avoid oscillations in OGS, we used finer meshes (about 10 K). The instabilities disappear and physical temperatures (between 0 and 1) can be obtained but the solution remains inaccurate (Fig. 8b). A careful check of the numerical parameters reveals that this inconsistency is related to the time stepping procedure. This procedure allows for using large time step sizes for which the solution becomes inaccurate.

To avoid this problem, OGS allows a maximum time step size. When this maximum step size is relatively small, OGS gives results which are in good agreement with the FG solution, as observed in Fig. 7. COMSOL allows for spurious oscillations to be removed by using a consistent stabilization technique that adds diffusion in the streamline direction. We used this technique with a mesh consisting of 10 K nodes, and we obtained a good agreement with the FG solution (Fig. 7). We should mention that the high dispersion cases are slightly sensitive to the computational mesh. However, for these cases, the computational efficiency of the solution becomes very dependent on the time integration scheme. In fact, for the high dispersion cases, the time required to reach the steady-state solution is much higher than the low dispersion cases. For longer transient simulations, we have observed that the higher-order integration schemes used in TRACES and COMSOL allow the use of large time step size and therefore reduce the computational cost. Table 1 also shows a good agreement between the non-dimensional metrics for Ra 1000 (errors are less than 3%). For Ra 5000 (low dispersion case), relatively weaker agreement can be observed between the different solutions, especially for the Nusselt number (errors are up to 13%). This could be attributed to the excessive numerical dispersion as for high Rayleigh numbers the convective term becomes dominant. The Nusselt numbers obtained with the three numerical codes are less than those obtained using the FG solution. This confirms that the discrepancy between the codes and FG solution is related to numerical dispersion. This numerical dispersion smooths out the temperature distribution and decreases the temperature gradient toward the hot wall.

The discrepancy between TRACES and FG is less than those for COMSOL and OGS. This Ra = 5,000 Ra =1,000

Low dispersion case

High dispersion case Ra =100 respectively. This is consistent with the results obtained in Fig. 6a which shows that the longitudinal dispersive flux is important in these zones as the velocity is parallel to the temperature gradient. A slight difference can be observed between the variance maps at Ra 1000 and 5000. This means that at high Rayleigh numbers the temperature distribution is mainly controlled by the dispersion coefficient. The behavior of the main isotherms as A L varies is depicted in Fig. 11a. This figure shows the main isotherms (0.1, 0.5 and 0.9) at Ra 100, R α disp 0.4 and for eight different values of A L ranging between 0.01 and 1. For the sake of brevity, we only present the results for Ra 100. Similar behaviors have been observed at Ra 1000 and 5000. Figure 11a shows that the 0.5 isotherm is more affected by A L than other isotherms. The increase of A L intensifies the mixing related to the thermal dispersion which widens the mixing zone. This pushes the hot isotherm (0.9) toward the hot wall at the top part of the cavity and leads to greater intrusion of isotherms in the lower part. Symmetrical behavior occurs for the cold isotherm (0.1). The 0.5 isotherm moves toward the hot wall (at the top) when A L is increased, as this latter dampens the convective flow. a more comprehensive understanding of thermal dispersion processes on NC. In the presence of thermal dispersion, the isotherms are different from those observed with the purely diffusive model. Their zigzag shape is mainly related to the velocity-dispersion processes.

The longitudinal dispersive flux is important toward the adiabatic walls where the velocity is parallel to the temperature gradient. Conversely, the longitudinal flux is important in front of the isothermal walls where the velocity is perpendicular to the temperature gradient.

A total of 144 runs of the FG method have been performed by varying the dispersion coefficients for three Rayleigh regimes (100, 1000 and 5000) in order to understand the effects of these dispersion parameters on fluid flow and heat transfer. In general, increasing thermal dispersion and we highlight the worthiness of these cases as benchmarks. We provide a new unambiguous benchmark that is more representative of physical processes than the purely diffusive one. From a numerical point of view, this work represents a generalization of the FG method to deal with natural convection in porous enclosures in the case of velocitydependent dispersion tensor. Furthermore, an important technical feature of this work is that the problem of natural convection in a square porous cavity can be solved in the spectral space using only temperature as primary unknown. We expect that this approach can be applied to other porous cavity geometries. The developed FG solution is limited to steady state. Extension of this work to transient cases is worthy of future inquiry.

The source code of the Fourier series solution and TRACES with several test examples are available on the author's Web site (https://lhyges.unistra.fr/FAHS-Marwan). The developed COMSOL and OGS models are uploaded on the software exchange Web sites.

Appendix A: Coefficients of the Spectral System

The coefficients of Eqs. ( 27)and( 28) are defined as follows: Gravity acceleration (m s -2 ) g 9.8

Viscosity (kg m -1 s -1 ) µ 10 -3

Fluid thermal conductivity (w m -1 °C-1 )0 . 6 5

Solid grains thermal conductivity (w m -1 °C-1 )1 . 5 9

Fluid specific heat (J kg °C-1 ) 4200

Solid grains specific heat (J kg °C-1 ) 580

Solid grains density 2000

Longitudinal dispersion coefficient (m) 1, 0 (for the high dispersion case) 0.001 (for the low dispersion case)

Transverse dispersion coefficient (m) 0, 1 (for the high dispersion case) 0.0001 (for the low dispersion case) See Table 3.