Noura Fajraoui

Marwan Fahs

Anis Younes

Bruno Sudret

Analyzing natural convection in porous enclosure with polynomial chaos expansions: Effect of thermal dispersion, anisotropic permeability and heterogeneity

Keywords:

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Natural convection (NC) in porous media can take place over a large range of scales that may go from fraction of centimeters in fuel cells to several kilometers in geological strata [START_REF] Nield | Convection in Porous Media[END_REF]. This phenomenon is related to the dependence of the saturating fluid density on the temperature 1 and/or compositional variations. A comprehensive bibliography about natural convection due to thermal causes can be found in the textbooks and handbooks by [START_REF] Nield | Convection in Porous Media[END_REF], [START_REF] Ingham | Transport phenomena in porous media III[END_REF], [START_REF] Vafai | Handbook of Porous Media[END_REF] and Vadasz [START_REF] Vadász | Emerging topics in heat and mass transfer in porous media: from bioengineering and microelectronics to nanotechnology[END_REF]. Comprehensive reviews on NC due to compositional effects have been provided by [START_REF] Diersch | Variable-density flow and transport in porous media: approaches and challenges[END_REF], [START_REF] Simmons | Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges[END_REF], [START_REF] Simmons | Variable density groundwater flow: From current challenges to future possibilities[END_REF] and Simmons et al [START_REF] Simmons | Variable density groundwater flow: From modelling to applications[END_REF]. NC in porous media can be encountered in a multitude of technological and industrial applications such as building thermal insulation, heating and cooling processes in solid oxide fuel cells, fibrous insulation, grain storage, nuclear energy systems, catalytic reactors, solar power collectors, regenerative heat exchangers, thermal energy storage, among others [START_REF] Nield | Convection in Porous Media[END_REF]; [START_REF] Ingham | Transport phenomena in porous media III[END_REF]; [START_REF] Ingham | Emerging Technologies and Techniques in Porous Media[END_REF]. Important applications can be also found in hydro-geology and environmental fields such as in geothermal energy Al-Khoury (2011); [START_REF] Carotenuto | A new methodology for numerical simulation of geothermal down-hole heat exchangers[END_REF], enhanced recovery of petroleum reservoirs [START_REF] Almeida | Integral transform methodology for convection-diffusion problems in petroleum reservoir engineering[END_REF]; [START_REF] Riley | Compositional variation in hydrocarbon reservoirs with natural convection and diffusion[END_REF]; [START_REF] Chen | Reservoir simulation: mathematical techniques in oil recovery[END_REF], geologic carbon sequestration [START_REF] Farajzadeh | Numerical simulation of density-driven natural convection in porous media with application for CO 2 injection projects[END_REF]; [START_REF] Class | A benchmark study on problems related to CO 2 storage in geologic formations[END_REF]; [START_REF] Islam | Numerical investigation of double diffusive natural convection of CO 2 in a brine saturated geothermal reservoir[END_REF][START_REF] Islam | Effects of geochemical reaction on double diffusive natural convection of CO 2 in brine saturated geothermal reservoir[END_REF]; [START_REF] Vilarrasa | Geologic carbon storage is unlikely to trigger large earthquakes and reactivate faults through which CO 2 could leak[END_REF], saltwater intrusion in coastal aquifers [START_REF] Werner | Seawater intrusion processes, investigation and management: recent advances and future challenges[END_REF] and infiltration of dense leachate from underground waste disposal [START_REF] Zhang | Multispecies contaminant plumes in variable density flow systems[END_REF].

Numerical simulation has emerged as a key approach to tackle the aforementioned applications in the last two decades. This is today a powerful and irreplaceable tool for understanding and predicting the behavior of complex physical systems. The literature concerning the numerical modeling and simulation of convective flow in porous media is abundant [START_REF] Holzbecher | Modeling density-driven flow in porous media: principles, numerics, software[END_REF]; [START_REF] Pop | Convective heat transfer: mathematical and computational modelling of viscous fluids and porous media[END_REF]; [START_REF] Viera | Mathematical and numerical modeling in porous media: applications in geosciences[END_REF]; [START_REF] Miller | Numerical simulation of water resources problems: Models, methods, and trends[END_REF]; [START_REF] Su | Modeling approaches to natural convection in porous media[END_REF]; [START_REF] Kolditz | Thermo-hydro-mechanical-chemical processes in fractured porous media: modelling and benchmarking[END_REF]. The NC in porous media is usually described by the conservation equations of fluid mass, linear momentum and energy, respectively. Either Darcy or Brinkman models are used as linear momentum conservation laws. Darcy model is a simplification of the Brinkman model by neglecting the effect of viscosity. This simplification is valid for low permeable porous media. For high permeable porous media Brinkman model is more suitable because the effective viscosity is about 10 times the fluid viscosity [START_REF] Givler | A determination of the effective viscosity for the Brinkman-Forchheimer flow model[END_REF]; [START_REF] Falsaperla | Double diffusion in rotating porous media under general boundary conditions[END_REF]; Shao et al. (2016). In the traditional modeling analysis of NC in porous media, the governing equations are solved under the assumption that all the parameters are known. However, in real applications, the determination of the input parameters may be difficult or inaccurate. For instance, in the simulation of geothermal reservoirs, the physical parameters (i.e. hydraulic conductivity and porosity) are subject to significant uncertainty because they are usually obtained by model calibration procedures, that are often carried out with relatively insufficient historical data O' Sullivan et al. (2001).

The uncertainties affecting the model inputs may have major effects on the model outputs. Typical examples about the significance of these effects (that are not exhaustive) can be found in the design of clinical devices or biomedical applications where small overheating can lead to unexpected serious disasters [START_REF] Davies | Sensitivity analysis of one-dimensional heat transfer in tissue with temperature-dependent perfusion[END_REF]; [START_REF] Ooi | Effects of natural convection within the anterior chamber on the ocular heat transfer[END_REF]; [START_REF] Wessapan | Aqueous humor natural convection of the human eye induced by electromagnetic fields: In the supine position[END_REF]. Hence, the evaluation of how the uncertainty in the model inputs propagates and leads to uncertainties in the model outputs is an essential issue in numerical modeling. In this context, uncertainty quantification (UQ) has become a must in all branches of science and engineering [START_REF] Brown | Assessing uncertainty propagation through physically based models of soil water flow and solute transport[END_REF]; [START_REF] Sudret | Uncertainty propagation and sensitivity analysis in mechanical models -Contributions to structural reliability and stochastic spectral methods[END_REF]; De [START_REF] De Rocquigny | Modelling under risk and uncertainty: an introduction to statistical, phenomenological and computational methods[END_REF]. It provides a rigorous framework for dealing with the parametric uncertainties. In addition, one wants to quantify how the uncertainty in the model outputs is due to the variance of each model input. This kind of studies is usually known as sensitivity analysis (SA) [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF]. UQ aims at quantifying the variability of a given response of interest as a function of uncertain input parameters, whereas GSA allows to determine the key parameters responsible for this variability. UQ and GSA are usually conducted by a multi-step analysis. The first step consists on the identification of model inputs that are uncertain and modelling them in a probabilistic context by means of statistical methods using data from experiments, legacy data or expert judgment. The second step consists in propagating the uncertainty in the input through the model. Finally, sensitivity analysis is carried out by ranking the input parameters according to their impact onto the prediction uncertainty. UQ and GSA have proven to be a powerful approach to assess the applicability of a model, for fully understanding the complex processes, designing, risk assessment and making decisions.

They have been extensively investigated in the literature (e.g. [START_REF] Saltelli | A quantitative model-independent method for global sensitivity analysis of model output[END_REF]; [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]; [START_REF] Xiu | A new stochastic approach to transient heat conduction modeling with uncertainty[END_REF]; [START_REF] Fesanghary | Design optimization of shell and tube heat exchangers using global sensitivity analysis and harmony search algorithm[END_REF]; [START_REF] Blackwell | A technique for uncertainty analysis for inverse heat conduction problems[END_REF]; [START_REF] Ghommem | Release of stored thermochemical energy from dehydrated salts[END_REF]; [START_REF] Sarangi | Manifold microchannel heat sink design using optimization under uncertainty[END_REF]; [START_REF] Zhao | Effect of parameters correlation on uncertainty and sensitivity in dynamic thermal analysis of thermal protection blanket in service[END_REF]; [START_REF] Mamourian | Two phase simulation and sensitivity analysis of effective parameters on turbulent combined heat transfer and pressure drop in a solar heat exchanger filled with nanofluid by response surface methodology[END_REF]; [START_REF] Shirvan | Numerical investigation and sensitivity analysis of effective parameters on combined heat transfer performance in a porous solar cavity receiver by response surface methodology[END_REF]; [START_REF] Rajabi | Sampling efficiency in Monte Carlo based uncertainty propagation strategies: application in seawater intrusion simulations[END_REF]). In the frame of flow and mass transfer in porous media, UQ and GSA have been applied to problems dealing with saturated/unsaturated flow [START_REF] Younes | Use of global sensitivity analysis to help assess unsaturated soil hydraulic parameters[END_REF], solute transport [START_REF] Fajraoui | Use of global sensitivity analysis and polynomial chaos expansion for interpretation of nonreactive transport experiments in laboratory-scale porous media[END_REF]; [START_REF] Ciriello | Polynomial chaos expansion for global sensitivity analysis applied to a model of radionuclide migration in a randomly heterogeneous aquifer[END_REF] and density driven flow Rajabi and Ataie-Ashtiani (2014); [START_REF] Riva | Probabilistic assessment of seawater intrusion under multiple sources of uncertainty[END_REF].

A careful literature review shows that the investigation of sensitivity analysis for NC in porous media has been limited to some special applications [START_REF] Shirvan | Numerical investigation and sensitivity analysis of effective parameters on combined heat transfer performance in a porous solar cavity receiver by response surface methodology[END_REF]. To the best of our knowledge, these analyses have never been performed for a problem involving NC within a porous enclosure. Yet, NC in porous enclosure has been largely investigated for several purposes [START_REF] Oztop | Investigation of natural convection in triangular enclosure filled with porous medi saturated with water near 4 •[END_REF]; [START_REF] Das | Studies on natural convection within enclosures of various (non-square) shapes-A review[END_REF] and several authors have contributed important results for such a configuration [START_REF] Bejan | On the boundary layer regime in a vertical enclosure filled with a porous medium[END_REF]; [START_REF] Prasad | Convective heat transfer in a rectangular porous cavity-effect of aspect ratio on flow structure and heat transfer[END_REF]; [START_REF] Beckermann | A numerical study of non-Darcian natural convection in a vertical enclosure filled with a porous medium[END_REF]; [START_REF] Gross | The application of flux-corrected transport (fct) to high Rayleigh number natural convection in a porous medium[END_REF]; [START_REF] Moya | Numerical study of natural convection in a tilted rectangular porous material[END_REF]; [START_REF] Lai | Natural convection across a vertical layered porous cavity[END_REF]; [START_REF] Baytaş | Entropy generation for natural convection in an inclined porous cavity[END_REF]; [START_REF] Saeid | Transient free convection in a square cavity filled with a porous medium[END_REF]; [START_REF] Saeid | Conjugate natural convection in a porous enclosure: effect of conduction in one of the vertical walls[END_REF]; [START_REF] Oztop | Investigation of natural convection in triangular enclosure filled with porous medi saturated with water near 4 •[END_REF]; [START_REF] Sojoudi | Unsteady natural convection within a porous enclosure of sinusoidal corrugated side walls[END_REF]; [START_REF] Chou | Effects of temperaturedependent viscosity on natural convection in porous media[END_REF]; [START_REF] Mansour | A numerical study on natural convection in porous media-filled an inclined triangular enclosure with heat sources using nanofluid in the presence of heat generation effect[END_REF].

Hence, keeping in view the various applications of NC in porous enclosure and the importance of uncertainty analysis in numerical modeling, a complete analysis involving GSA and UQ study is developed in this work to address this gap. The considered problem deals with the square porous cavity. Such a problem is widely used as a benchmark for numerical code validation due to the simplicity of the boundary conditions [START_REF] Walker | Convection in a porous cavity[END_REF]; [START_REF] Manole | Numerical benchmark results for natural convection in a porous medium cavity[END_REF]; [START_REF] Misra | A comparative study of porous media models in a differentially heated square cavity using a finite element method[END_REF]; [START_REF] Baytaş | Entropy generation for natural convection in an inclined porous cavity[END_REF]; [START_REF] Alves | Transient natural convection inside porous cavities: hybrid numerical-analytical solution and mixed symbolic-numerical computation[END_REF]; Fahs et al. (2015a); Shao et al. (2016,?); [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 physical insights and better understanding of NC processes in porous media [START_REF] Getachew | Natural convection in a porous cavity saturated with a non-newtonian fluid[END_REF]; [START_REF] Baytaş | Entropy generation for natural convection in an inclined porous cavity[END_REF]; [START_REF] Saeid | Transient free convection in a square cavity filled with a porous medium[END_REF]; [START_REF] Leong | Natural convection in rectangular layered porous cavities[END_REF]; [START_REF] Mahmud | Mixed convection in a square vented enclosure filled with a porous medium[END_REF]; [START_REF] Choukairy | Numerical and analytical analysis of the thermosolutal convection in an heterogeneous porous cavity[END_REF]; [START_REF] Malashetty | Linear and nonlinear double-diffusive convection in a fluidsaturated porous layer with cross-diffusion effects[END_REF] . As model inputs, we consider the physical parameters characterizing the porous media and the saturating fluid as the permeability, porosity, thermal diffusivity and thermal expansion. All these parameters can be described by the Rayleigh number which represents the ratio between the buoyancy and diffusion effects.

A common simplification for NC in porous media is to consider the saturated porous media with an equivalent thermal diffusivity (based on the porosity) and neglect the key process of heat mixing related to velocity dependent dispersion. Yet several studies have found that thermal dispersion plays an important role in NC systems [START_REF] Howle | Natural convection in porous media with anisotropic dispersive thermal conductivity[END_REF]; [START_REF] Metzger | Optimal experimental estimation of thermal dispersion coefficients in porous media[END_REF]; [START_REF] Pedras | Thermal dispersion in porous media as a function of the solid-fluid conductivity ratio[END_REF]; [START_REF] Jeong | Estimation of the thermal dispersion in a porous medium of complex structures using a lattice boltzmann method[END_REF]; [START_REF] Yang | Analysis of heat flux bifurcation inside porous media incorporating inertial and dispersion effects-an exact solution[END_REF]; [START_REF] Kumar | Natural convection in an anisotropic porous enclosure due to nonuniform heating from the bottom wall[END_REF]; [START_REF] Sheremet | Effect of thermal dispersion on transient natural convection in a wavy-walled porous cavity filled with a nanofluid: Tiwari and das nanofluid model[END_REF]; [START_REF] Plumb | Non-Darcy natural convection from heated surfaces in saturated porous media[END_REF]; [START_REF] Cheng | Thermal dispersion effects in non-Darcian convective flows in a saturated porous medium[END_REF]; [START_REF] Hong | Analysis of thermal dispersion effect on vertical-plate natural convection in porous media[END_REF]; [START_REF] Cheng | Transverse thermal dispersion and wall channelling in a packed bed with forced convective flow[END_REF]; [START_REF] Amiri | Analysis of dispersion effects and non-thermal equilibrium, non-Darcian, variable porosity incompressible flow through porous media[END_REF]; [START_REF] Shih-Wen | Non-uniform porosity and thermal dispersion effects on natural convection about a heated horizontal cylinder in an enclosed porous medium[END_REF] and applications related to transport in natural porous media [START_REF] Abarca | Anisotropic dispersive Henry problem[END_REF]; [START_REF] Jamshidzadeh | Fluid dispersion effects on density-driven thermohaline flow and transport in porous media[END_REF]; [START_REF] Fahs | The Henry problem: new semi-analytical solution for velocity-dependent dispersion[END_REF]. Hence, main attention is given here to understand the impact of anisotropic thermal dispersion by including the longitudinal and transverse dispersion coefficients in the model inputs. Furthermore, anisotropy in the hydraulic conductivity is acknowledged as it is one of the properties of porous media which is a consequence of asymmetric geometry and preferential orientation of the solid grains Shao et al. (2016).

Finally heterogeneity of the porous media is considered as a source of uncertainty as it has a significant impact on NC in porous media [START_REF] Simmons | Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges[END_REF]; [START_REF] Nield | A discussion on the effect of heterogeneity on the onset of convection in a porous medium[END_REF]; [START_REF] Kuznetsov | Natural convective boundary-layer flow of a nanofluid past a vertical plate[END_REF]; [START_REF] Zhu | Entropy generation due to three-dimensional double-diffusive convection of power-law fluids in heterogeneous porous media[END_REF]. As model outputs, we consider different quantities that are often used to assess the flow and the heat transfer processes in porous cavity as the temperature spatial distribution, the Nusselt number and the maximum velocity components.

In this work, we perform a global sensitivity analysis using a variance-based technique. [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models[END_REF]; [START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF]. Monte Carlo (MC) is one of the most commonly used methods. However, it might become impractical, because of the large number of repeated simulations required to attain statistical convergence of the solution, especially for complex problem (e.g., [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]; [START_REF] Ballio | Convergence assessment of numerical monte carlo simulations in groundwater hydrology[END_REF] and references therein). In this context, new approaches based on advanced sampling strategies have been introduced to reduce the computational burden associated with Monte Carlo simulations. Among different alternatives, Polynomial Chaos Expansions (PCE) have been shown to be an efficient method for UQ and GSA Blatman andSudret (2010b,a, 2011). In PCE, the key idea is to expand the model response in terms of a set of orthonormal multivariate polynomials orthogonal with respect to a suitable probability measure [START_REF] Ghanem | Stochastic finite elements -A spectral approach[END_REF]. They allow one to uncover the relationship between different input parameters and how they affect the model outputs. Once a PCE representation is available, the Sobol' sensitivity indices are then obtained via a straightforward post-processing analysis without any additional computational cost [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. It can also be used to perform an uncertainty quantification using Monte Carlo analysis at a significantly reduced computational cost (see, e.g., [START_REF] Fajraoui | Use of global sensitivity analysis and polynomial chaos expansion for interpretation of nonreactive transport experiments in laboratory-scale porous media[END_REF] and references therein).

The structure of the present study is as follows. Section 2 is devoted to the description of the benchmark problem and the governing equations. Section 3 describes the numerical model. Section 4 describes the sensitivity analysis procedure using Sparse PCE. Section 5 discusses the GSA and UQ results for homogeneous and heterogeneous porous media. Finally, a summary and conclusions are given in Section 6.

Problem statement and mathematical model

The system under consideration is a square porous enclosure of length H filled by a saturated heterogeneous porous medium. The properties of the fluid and the porous medium are assumed to be independent on the temperature. The porous medium and the saturating fluid are locally in thermal equilibrium. We assume that the Darcy and Boussinesq approximations are valid and that the inertia and the viscous drag effects are negligible. Under these conditions, the fluid flow in anisotropic porous media can be described using the continuity equation and Darcy's law written in Cartesian coordinates as follows [START_REF] Mahmud | Mixed convection in a square vented enclosure filled with a porous medium[END_REF]:

∂u ∂x + ∂v ∂y = 0, (1)
u = - k x µ ∂p ∂x , (2) v = - k y µ ∂p ∂y + (ρ -ρ c ) g , (3)
where u and v LT -1 are the fluid velocity components in the x and y directions; p M L -1 T -2 is the total pressure (fluid pressure and gravitational head); k x and k y L 2 are the permeability components in the x and y directions; µ M L -1 T -1 is the dynamic viscosity; ρ and ρ c M L -1 being respectively the density of the mixed fluid and density of the cold fluid; and g LT -2 is the gravitational constant.

The heat transfer inside the cavity is modeled using the energy equation written as:

∂T ∂t +u ∂T ∂x +v ∂T ∂y = α m ∂ 2 T ∂x 2 + ∂ 2 T ∂y 2 + ∂ ∂x α xx disp ∂T ∂x + α xy disp ∂T ∂y + ∂ ∂y α xy disp ∂T ∂x + α yy disp ∂T ∂y . (4)
Here T [Θ] is the temperature; α m L 2 T -1 is the the effective thermal diffusivity and

α disp L 2 T -1
is the thermal dispersion tensor. In this work, we use the nonlinear model with anisotropic tensor as in [START_REF] Howle | Natural convection in porous media with anisotropic dispersive thermal conductivity[END_REF], that is defined as follows:

α xx disp = (α L -α T ) u 2 √ u 2 + v 2 + α T u 2 + v 2 , (5)
α xy disp = (α L -α T ) u • v √ u 2 + v 2 , (6)
α yy disp = (α L -α T ) v 2 √ u 2 + v 2 + α T u 2 + v 2 , ( 7
)
where α L and α T [L] are respectively the longitudinal and transversal dispersivity coefficient, which are considered uniform in the system. No-flow boundary conditions are assumed across all boundaries. The left and right vertical walls are maintained at constant temperatures T h and T c such that (T h > T c ), respectively. The horizontal surfaces are assumed to be adiabatic (Fig. 1).

The system (1)-( 7) is completed by specifying a constitutive relationship between fluid properties ρ and the temperature T . The density of the mixed fluid is assumed to vary with temperature as a first-order polynomial, that is:

ρ = ρ c [1 -β(T -T c )], ( 8
)
where β is the coefficient of thermal expansion.

Special attention is given in the literature to the NC in heterogeneous porous media because of the nonuniformity of the permeability and/or the thermal diffusivity affect significantly the overall rate of the heat transfer. The effect of heterogeneity is especially challenging in geothermal applications since hydraulic properties, such as permeability, can vary by several orders of magnitude over small spatial scales Fahs et al. (2015a). In this context, the impact of heterogeneity has been studied for both external and internal natural convection using different heterogeneity configurations (stratified, horizontal, vertical, random, and periodic) [START_REF] Kuznetsov | Natural convective boundary-layer flow of a nanofluid past a vertical plate[END_REF] and references therein. In this work, the heterogeneity of the porous media is described via the exponential model as in Fahs et al. x-and y-directions are given by:

k x = k x0 e σy , (9)
k y = k y0 e σy , (10)
where k x0 and k y0 are respectively the permeability at y = 0 and x = 0; σ is the rate of change of ln(k) in the y direction.

The system (1)-(3) can be reformulated in dimensionless form as:

∂u * ∂x * + ∂v * ∂y * = 0, (11)
u * = -e σ * y * ∂p * ∂x * , (12)
v * = -r k e σ * y * ∂p * ∂y * + Ra • T * , ( 13
)
where

p * = k x0 µα m p; u * = H α m u; v * = H α m v; T * = T -T c ∆T ; ∆T = T h -T c
is the temperature difference between hot and cold walls;

x * = x H ; y * = y H ;σ * = σ • H; r k = k y0 k x0 . Ra represents the
Rayleigh number, that is given by:

Ra = k y • ρ c • β • g • ∆T • H µα m . ( 14
)
The steady-state energy equation can be written in dimensionless form as:

u * ∂T * ∂x * +v * ∂T * ∂y * = ∂ 2 T * ∂x * 2 + ∂ 2 T * ∂y * 2 + ∂ ∂x * α xx,
α xx, * disp = (α * L -α * T ) (u * ) 2 (u * ) 2 + (v * ) 2 + α T (u * ) 2 + (v * ) 2 , ( 16
)
α xy, * disp = (α * L -α * T ) u * • v * (u * ) 2 + (v * ) 2 , ( 17
)
α yy, * disp = (α * L -α * T ) (v * ) 2 (u * ) 2 + (v * ) 2 + α * T (u * ) 2 + (v * ) 2 , ( 18
)
where

α * L = α L H ; α * T = α T H .

The numerical model

Numerical simulation of thermal-driven transfer problem is highly sensitive to discretization errors. Furthermore, hydraulic anisotropy, heterogeneity and anisotropic thermal dispersion render the numerical solution more challenging as they require specific numerical techniques. Therefore, it is extremely important to select the appropriate numerical methods for solving the governing equations. In this work, we use the advanced numerical model developed by [START_REF] Younes | Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods[END_REF] and Younes and Ackerer [START_REF] Younes | Solving the advection-dispersion equation with discontinuous galerkin and multipoint flux approximation methods on unstructured meshes[END_REF]. In this model, appropriate techniques for both time integration and spatial discretization are used to simulate coupled flow and heat transfer. For the spatial discretization, a specific method is used to achieve high accuracy for each type of equation. Thus, the Mixed Hybrid Finite method is used for the discretization of the flow equation. This method produces accurate and consistent velocity fields even for highly heterogeneous domains [START_REF] Farthing | Mixed finite element methods and higher-order temporal approximations[END_REF]; [START_REF] Durlofsky | Accuracy of mixed and control volume finite element approximations to Darcy velocity and related quantities[END_REF]. The heat transfer equation is discretized through a combination of a discontinuous Galerkin (DG) and Multipoint flux approximation (MPFA) methods. For the convective part, the DG method is used because it provides robust and accurate numerical solutions for problems involving steep fronts [START_REF] Younes | Solving the advection-dispersion equation with discontinuous galerkin and multipoint flux approximation methods on unstructured meshes[END_REF]; [START_REF] Tu | A slope limiting procedure in discontinuous galerkin finite element method for gasdynamics applications[END_REF]. For the diffusive part, the MPFA method is used because it allows for the handling of anisotropic heterogeneous domains and can be easily combined with the DG method [START_REF] Younes | Solving the advection-dispersion equation with discontinuous galerkin and multipoint flux approximation methods on unstructured meshes[END_REF]. The method of lines (MOL) is used for the time integration. This method improves the accuracy of the solution through the use of adaptive higherorder time integration schemes with formal error controls. The numerical model has been validated against lab experimental data for variable density flow [START_REF] Konz | Variable-density flow in heterogeneous porous media: Laboratory experiments and numerical simulations[END_REF]. It has been also validated by comparison against semi-analytical solutions for NC in porous square cavity Fahs et al. (2014,?), seawater intrusion in heterogeneous coastal aquifer [START_REF] Younes | Extension of the Henry semi-analytical solution for saltwater intrusion in stratified domains[END_REF] and seawater intrusion in anisotropic dispersive porous media [START_REF] Fahs | The Henry problem: new semi-analytical solution for velocity-dependent dispersion[END_REF]. We should note that the numerical model allows for transient simulations while in the GSA study we consider the steady state solutions. Hence transient solutions are performed until a long nondimensional time to ensure steady conditions

Polynomial Chaos Expansion for Sensitivity Analysis

GSA is a useful tool that aims at quantifying which input parameters or combinations thereof contribute the most to the variability of the model responses, quantified in terms of its total variance. Variance-based sensitivity method have gained interest since the mid 90's, in this particular context.

Here, we base our analysis on the the Sobol' indices, which are widely used as sensitivity metrics Sobol' (1993) and do not rely on any assumption regarding the linearity or monotonous behavior of the physical model.

In the sequel, we consider Y = M (X), a mathematical model that describe a scalar output of the considered physical system, which depends on M -uncertain input parameters. M may, represent a scalar model response. In the case of vector-valued response, i.e.; {Y ∈ R N , N > 1}, the following approach may be applied component-wise. We consider the uncertain parameters as independent random variables gathered into a random vector X = {X 1 , ..., X M } with joint probability density function (PDF) f X and marginal PDFs {f Xi (x i ), i = 1, ..., M }. Within this context, we will introduce next the variance-based Sobol' indices. The interested reader is referred to [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]; [START_REF] Gratiet | Metamodel-based sensitivity analysis: Polynomial chaos expansions and Gaussian processes[END_REF] for a deeper insight into the details.

Anova-based sensitivity indices

Provided that the function M is square-integrable with respect to the probability measure associated with f X , it can be expanded in summands of increasing dimension as:

M(X) = M 0 + M i=1 M i (X i ) + 1≤i<j≤M M ij (X i , X j ) + ... + M 12...M (X), ( 19
)
where M 0 is the expected value of M(X), and the integrals of the summands M i1,i2,...,is with respect to their own variables is zero, that is:

DX i k M i1,i2,...,is (X i1,i2,...,is )f Xi k (X i k ) = 0 for 1 ≤ k ≤ s, (20)
where D Xi k and f Xi k (X i k ) respectively denote the support and marginal PDF of X i k .

Eq (19) can be written equivalently as:

M(X) = M 0 + u =0 M u (X u ). ( 21
)
Here u = {i 1 , i 2 , ..., i M } ⊆ {1, 2, ..., M } are index sets and X u are subvectors containing only those components of which the indices belong to u. This representation is called the Sobol' decomposition.

It is unique under the orthogonality conditions between summand, namely:

E[M u (X u )M v (X v )] = 0. ( 22
)
Thanks to the uniqueness and orthogonality properties, it is straightforward to decompose the total variance of Y , denoted D in a sum of partial variance D u :

D = Var [M(X)] = u =0 D u = u =0 Var [M u (X u )] , (23)
where:

D u = Var [M u (X u )] = E[M 2 u (X u )]. (24)
This leads to a natural definition of the Sobol' Indices S u :

S u = D u /D, (25)
which measures the amount of the total variance due to the contribution of the subset X u . In particular, the first-order sensitivity index is defined by:

S i = D i /D. ( 26
)
The first-order sensitivity indices S i measures the amount of variance of Y that is due to the parameter X i considered separately. The overall contribution of a parameter X i to the response variance including its interactions with the other parameters is then given by the total sensitivity indices. They include the main effects S i and all the joint terms involving parameter X i , i.e.

S T i = Ii D u /D, I i = {u ⊃ i}, (27)
In principle, one should rely upon the total sensitivity index to infer the relevance of the parameters [START_REF] Saltelli | On the relative importance of input factors in mathematical models: safety assessment for nuclear waste disposal[END_REF]. The higher S T i , the more X i is an important parameter for the model response. In contrast, X i is termed unimportant (in terms of probabilistic modelling) if S T i = 0.

The evaluation of Sobol' indices requires the computation of 2 M Monte Carlo integrals of the model response M (X). This can be costly to manage, especially when dealing with time-consuming computational models. Fortunately, the Sobol' indices can easily be computed using the Polynomial Chaos Expansion (PCE) technique [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. They are analytically obtained via a straightforward post-processing of the expansion. The PCE will be described in the next section.

Polynomial chaos expansion

The model response can be casted into a set of orthonormal multivariate polynomial as:

Y = M(X) = α∈A y α Ψ α (X), ( 28
)
where A is a multi-index α = {α 1 , ..., α M }, {y α , α ∈ A} are the expansion coefficients to be determined, {Ψ α (X), α ∈ A} are multivariate polynomials which are orthonormal with respect to the joint

pdf f X of X, i.e E[Ψ α (X)Ψ β (X)] = 1 if α = β and 0 otherwise.
The multivariate polynomials Ψ α are assembled as the tensor product of their appropriate univariate polynomials, i.e

Ψ α (x) = M i=1 φ (i) αi (x i ), (29)
where φ

(i)
αi is a polynomial in the i-th variable of degree α i . These bases are chosen according to the distributions associated with the input variables. For instance, if the input random variables are standard normal, a possible basis is the family of multivariate Hermite polynomials, which are orthogonal with respect to the Gaussian measure. Other common distributions can be used together with basis functions from the Askey scheme [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]. A more general case can be treated through an isoprobabilistic transformation of the input random vector X into a standard random vector. The set of multi-indices A in Eq. ( 28) is determined by an appropriate truncation scheme. In the present study, a hyperbolic truncation scheme [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] is employed, which consists in selecting all polynomials satisfying the following criterion:

|α q = M i=1 α q i 1/q ≤ p, (30)
with p being the highest total polynomial degree, 0 < q ≤ 1 being the parameter determining the hyperbolic truncation surface. This truncation scheme allows for retaining univariate polynomials of degree up to p, whereas limiting the interaction terms.

The next step is the computation of the polynomial chaos coefficients {y α , α ∈ A}. Several intrusive (e.g. Galerkin scheme) or non-intrusive approaches (e.g. stochastic collocation, projection, regression methods) [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]; [START_REF] Xiu | Numerical methods for stochastic computations: a spectral method approach[END_REF] are proposed in the literature. We herein focus our analysis on the regression methods also known as least-square approaches. A set of N realization of the input vector, X = {x (1) , ..., x (N ) }, is then needed, called experimental design (ED). The set of coefficient are then computed by means of the least-square minimization method, that is:

ŷα = argmin y α ∈R cardA 1 N N i=1 M x (i) - α∈A y α Ψ α (x (i) ) 2 . ( 31
)
The number of terms in Eq. ( 28) may be unnecessarily large, thus a sparse PCE can be more efficient to capture the behavior of the model by disregarding insignificant terms from the set of regressors. We herein adopt the least angle regression (LAR) method proposed in [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] which involves a sparse representation containing only a small number of regressors compared to the classical full representation. The reader is referred to [START_REF] Efron | Least angle regression[END_REF] for more details on the LARS technique and to [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on Least Angle Regression[END_REF] for its implementation in the context of adaptive sparse PCE.

It can be worth noting that the constructed PCE can also be employed as a surrogate model of the target output in cases when evaluating a large number of model responses is not affordable. It is thus important to assess its quality. A good measure of the accuracy is the Leave-One-Out (LOO) error, which allows a fair error estimation at an affordable computational cost Blatman and Sudret (2010a).

The relative LOO error is defined as:

ǫ LOO = N i=1 M(x (i) ) -M P C (x (i) ) 1 -h i 2 N i=1 M(x (i) ) -μY 2 , ( 32
)
where h i is the i th diagonal term of matrix

Ψ(Ψ T Ψ) -1 Ψ T , where Ψ = {Ψ ij = Ψ j X i } and μY = 1 N N i=1 M(x (i) ).

Polynomial chaos expansions for sensitivity analysis

Once the PCE is built, the mean µ and the total variance D can be obtained using properties of the orthogonal polynomials [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF], such that:

µ = y 0 , (33)
D = α∈A\0 y 2 α . ( 34
)
As mentioned above, the Sobol' indices of any order can be computed in a straightforward manner. The first order and total Sobol' indices are then given by [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]:

S i = α∈Ai y 2 α /D, A i = {α ∈ A : α i > 0, α j =i = 0}, (35)
and

S T i = α∈A T i y 2 α /D, A T i = {α ∈ A : α i > 0}. ( 36
)
Of particular interest is the marginal effect (also called univariate effect, see Deman et al. ( 2016))

of the parameters X i , which enables investigation of the range of variation across which the model response is most sensitive to X i . It corresponds to the sum of the mean values and first-order summands comprising univariate polynomials only, i.e.:

E[M(X) | X i = x i ] = M 0 + α∈Ai y α Ψ α (x i ). ( 37
)

Results and discussions

The PCEs presented in the previous section are used to perform UQ and GSA for the problem of natural convection in porous square cavity. The dimensionless form of the governing equations leads to define the model input parameters as follows:

• The average Rayleigh number (Ra): the Rayleigh number represents the ratio between the buoyancy and the diffusion effects. It depends on the porous media properties (porosity, thermal diffusivity and permeability), fluid properties (thermal diffusivity, viscosity, density and thermal expansion), the characteristic domain length and the temperature gradient. For isotropic porous media, the Rayleigh number is defined based on the scalar permeability of the porous media. For the general case of an anisotropic porous media, Ra is defined based on the permeability in the vertical direction (k y ) [START_REF] Bennacer | Double diffusive convection in a vertical enclosure filled with anisotropic porous media[END_REF]. In this work, we are concerned with anisotropic heterogeneous porous media. Thus, we distinguish between local Rayleigh number based on the local permeability (see Eq. ( 14)) and the average Rayleigh number based on the overall average permeability Fahs et al. (2015a). The local Ra number can be formulated as follows:

Ra = Ra 0 e σ * y * , ( 38
)
where Ra 0 is the local Rayleigh number at the bottom of the domain. The average Rayleigh number Ra is then obtained by integrating the local Ra overall the domain, that is given by:

Ra = Ra 0 1 0 e σ * y * dy * = e σ * -1 σ * Ra 0 . ( 39
)
The range of variability of the average Rayleigh number Ra is from 0 to 1000. This range of variation is physically plausible.

• The permeability anisotropy ratio (r k ): this ratio is commonly used to describe the hydraulic anisotropy of the porous media [START_REF] Abarca | Anisotropic dispersive Henry problem[END_REF]; [START_REF] Bennacer | Double diffusive convection in a vertical enclosure filled with anisotropic porous media[END_REF]. In this work, the model used to describe the heterogeneity of the porous media leads to a constant anisotropy ratio, calculated based on the permeability on the x and y directions at the bottom of the domain (k x0 and k y0 ). As a common practice in porous media, the range of variability of r k is considered to be between 0 and 1 [START_REF] Abarca | Anisotropic dispersive Henry problem[END_REF].

• The non-dimensional dispersion coefficients (α * L and α * T ): these parameters correspond to the longitudinal and transverse thermal dispersion coefficients. They account for the enhancement of heat transfer due to hydrodynamic dispersion. The longitudinal dispersion α * L corresponds to the heat transfer along the local (Darcy) velocity vector while the transverse dispersion α * T acts normally to the local velocity. A detailed review about the physical understanding of longitudinal and transverse thermal dispersion is given in Howle and Georgiadis [START_REF] Howle | Natural convection in porous media with anisotropic dispersive thermal conductivity[END_REF].

According to Howle and Georgiadis [START_REF] Howle | Natural convection in porous media with anisotropic dispersive thermal conductivity[END_REF], Abarce et al., [START_REF] Abarca | Anisotropic dispersive Henry problem[END_REF] and [START_REF] Fahs | The Henry problem: new semi-analytical solution for velocity-dependent dispersion[END_REF], α * L was varied between 0.1 and 1 and α * T between 0.01 and 0.1.

• The rate of heterogeneity variation (σ * z ): this parameter is used to quantify the effect of the heterogeneity distribution on the model outputs. In fact, the geometrical distribution of the heterogeneity in the computational domain is not often well-defined. For instance, in hydrogeology, the heterogeneity distribution cannot be clearly described because hydraulic parameters do not correlate well with lithology. As in Fahs et al., Fahs et al. (2015a) and shao et al., Shao et al. (2016), the range of variability of (σ * z ) is assumed to be from 0 to 4. Uncertainty in these parameters is related to our imperfect knowledge of the porous media properties (porosity, permeability tensor, heterogeneity distribution) and thermo-physical parameters of both porous media grains and saturating fluid (thermal conductivity and dispersion). Without further information, and in view of drawing general conclusions, uniform distributions are selected for all parameters. Moreover, the parameters are assumed to be statistically independent.

The results of the numerical model will be analyzed using several quantities of interest (QoI) which are controlled by the model inputs. To describe flow process we use the maximum dimensionless velocity components (u * max and v * max ). For the heat transfer process, the assessment is based on the spatial distribution of the dimensionless temperature (T * ). In addition, and as it is customary for the cavity problem, the heat processes are assessed using the wall average Nusselt number N u given by:

N u = 1 0 N u(y * ) dy * , ( 40
)
where N u is the local Nusselt number. The local Nusselt number represents the net dimensionless heat transfer at a local point on the hot wall. It is defined as the ratio of the total convective heat flux to its value in the absence of convection. When thermal dispersion is considered, the local Nusselt number is defined as follows [START_REF] Howle | Natural convection in porous media with anisotropic dispersive thermal conductivity[END_REF]; [START_REF] Sheremet | Effect of thermal dispersion on transient natural convection in a wavy-walled porous cavity filled with a nanofluid: Tiwari and das nanofluid model[END_REF]:

N u = 1 + α * T (u * ) 2 + (v * ) 2 ∂T * ∂x * x * =0 . ( 41
)
5.1 Homogeneous case

Numerical details

Preliminary simulations were performed for different grid size in order to test the influence of grid discretization on the QoIs. They were performed under regular triangular mesh obtained by subdividing square elements into four equal triangles (by connecting the center of each square to its four nodes).

Regular grids are used here to avoid instabilities and inaccuracies that can be caused by the change in mesh sizes within irregular grids. The most challenging configuration of the uncertain parameters is considered. This corresponds to the case with the highest Rayleigh number (Ra = 1, 000) and anisotropy ratio (r k = 1) and lowest values of longitudinal and transverse thermal dispersion coefficients. For such a case, the heat transfer process is mainly dominated by the buoyancy effects which are at the origin of the rotating flow within the cavity. As a consequence, the steady state isotherms are sharply distributed and they have a spiral shape as they follow the flow structure due to the small thermal diffusivity. A relatively fine mesh should be used in this case to obtain a mesh independent solution. Several simulations are performed by increasing progressively the mesh refinement and by comparing the solution for two consecutive levels of grid. The tests revealed that the uniform grid formed by 40, 000 elements is adequate to render accurate results and capture adequately the flow and heat transfer processes. All simulations were run for 8 minutes because this is the required time for the homogeneous problem to reach the steady state solution. These discretization parameters are kept fixed in subsequent simulations.

In view of computing the PCE expansion of the model outputs in terms of the 4 input random variables X = {Ra, r k , α * L , α * T }, several sets of parameter values sampled according to their respective pdf's are needed. For this purpose, we use an experimental design of size N = 150 drawn with Quasi Monte Carlo sampling (QMC). It is a well-known technique for obtaining deterministic experimental designs that covers at best the input space ensuring uniformity of each sample on the margin input variables. In particular, Sobol' sequences are used. PCE meta-models is constructed by applying the procedure described in Section 3 for the considered QoI's. In the case of multivariate output (temperature), a PCE is constructed component-wise (i.e.; for each points of the grid). The candidate basis is determined using a standard truncation scheme (see Eq. ( 27)) with q = 1. The maximum degree p is varied from 1 to 20 and the optimal sparse PCE is selected by means of the corrected relative LOO error (see Eq. ( 29)). The corresponding results (e.g. polynomial degree giving the best accuracy, relative LOO error and number of retained polynomials) of the PCE are given in Table 1 for the three scalar output N u, u * max and v * max . For instance, when the average Nusselt number is considered, the optimal PCE is obtained for p = 4 and the corresponding LOO error is err LOO = 8.6 × 10 -4 . The sparse meta-model includes 63 basis elements, whereas the size of full basis is 70. In Fig. 2, we compare the values of the PCE with the respective values of the physical model at a validation set consisting of 1, 000 MC simulations. We note that these simulations do not coincide with the ED used for the construction of the PCE. An excellent match is observed for both N u and v * max ; which is also illustrated by a small LOO error (less than 0.001). Discrepancies between PCE and true model are observed for u * max ; especially for larger values of u * max . The related LOO error is equal to 5.81 × 10 -2 .

Global sensitivity analysis

This section is devoted to GSA in order to identify the most influential parameters and to understand the marginal effect of the parameters onto the model outputs. Depending on the output QoI's, a 1, the relative LOO error varies from 0.04% to 5.8%. It is important to emphasize that excellent GSA results are obtained by PCE as soon as err LOO < 10 -3 . Moreover, the results obtained for u * max are also deemed acceptable.

• GSA of the temperature distribution Fig. 3a illustrates the spatial distribution of the mean of the temperature based on PCE. In this case, the presented approach is applied component-wise. Indeed, the PCE of a numerical model with many outputs is carried-out by metamodelling independently each model output. Fig. 3a shows that the distribution of the mean temperature reflects the general behavior of the heat transfer in the case of NC in square porous cavity. The isotherms are not vertical as they are affected by the circulation of the fluid saturating the porous media. In order to evaluate how far the temperatures are spread out from their mean, we plot in Fig. 3b L and α * T . We recall that the total Sobol' indices involve the total effect of a parameter including nonlinearities as well as interactions. Thus, they allow us to rank the parameters according to their importance. Focusing on Fig. 4, we can see that the most influential parameters are α * L and α * T . A complementary effect between these parameters is observed. The effect of α * L is more pronounced than that for α * T as its zone of influence is located in the region where the temperature variance is maximum. It is worth nothing that complementary effect between the influence of α * L and α * T can be explained by reformulating the dispersive heat flux in terms of the dot product of the velocity and temperature gradient vectors. Considering Eqs. ( 5)-( 7) the thermal dispersive flux q disp can be rearranged as follows:

q x disp q z disp =     α T | V | ∂T ∂x - u | V | (V • ∇T ) + uα L | V | (V • ∇T ) α T | V | ∂T ∂y - v | V | (V • ∇T ) + vα L | V | (V • ∇T )     . ( 42
)
where V is the velocity vector. This equation reveals that the longitudinal (resp. transverse) dispersion is an increasing (resp. decreasing) function of V • ∇T . Around the top and bottom surfaces of the cavity, the velocity is almost horizontal and parallel to the thermal gradient. Hence, V • ∇T exhibits its maximum value and by consequence temperature distribution in these zones is mainly controlled by α * L (see Fig. 4). The velocity near the vertical walls is vertical and relatively perpendicular to the thermal gradient. V • ∇T tends towards zero (minimum value). Hence, the temperature distribution in these regions is mainly controlled by α * T (see Fig. 4d).

The sensitivity of the temperature field due to Ra is less important than that of α * L and α * T (Fig. 4a). Its zone of influence expands along the cavity diagonal bisector with increasing magnitude toward the cavity center and near the corners. The effect of Ra on the temperature distribution is related to the heat mixing by thermal diffusion and/or the convection due to the buoyancy effects. Around the top right and bottom left corners the fluid velocity is very small as it should be zero right on the corners to satisfy the boundary conditions. Hence heat mixing between hot and cold fluids by thermal dispersion is almost negligible. In addition, around these corners the thermal boundary layers are very thin. Thus, the temperature gradient is very important so that mixing by thermal diffusion is dominating and by consequence temperature distribution is sensitive to Ra. In the center of the cavity the thermal dispersion tensor is also negligible because the rotating flow is relatively slow. Hence mixing is mainly related to diffusion and by consequence sensitivity to the Ra is relatively important. Fig. 4b indicates that the temperature distribution is slightly sensitive to the permeability ratio. The zone of influence of r k matches well with the region in which the flow is strongly bidirectional. This is physically understandable, since r k expresses the ability of a porous media to transmit fluid in a direction perpendicular to the main flow. Hence r k is a non-influent parameter in the zones where the flow is almost unidirectional.

• GSA of the scalar QoIs Fig. 5 shows bar-plots of the first order and total Sobol indices of the average Nusselt number N u. Inspection of the sensitivity indices showed that the variability of N u is mainly due to the principal effects of Ra and α * T . The most influential parameter is Ra with S T Ra = 0.8. A small influence of r k and α * L is also observed. Interactions between the random parameters are not significant (not shown). The maximum value obtained is S Ra,α * T = 0.035. The results are summarized in Table 2. To further elaborate our investigation, we examine the marginal effect of the uncertain parameters on the model response. This effect corresponds to the evolution of the model output with respect to a single parameter averaged on the other parameters (Eq. ( 37)). Indeed, if a parameter is sensitive, Fig. 6a demonstrates that N u increases with Ra. Indeed, the increase of Ra enhances the buoyancy effects and reduces the thickness of the thermal boundary layer in the lower part of the hot wall. This leads to higher values of N u as the temperature distribution becomes steep near the hot wall, especially around the bottom corner. Same behavior of N u is observed against α * T (Fig. 6d). This parameter affects slightly the velocity field (as it will be shown later in this paper). It slightly affects the temperature distribution at the hot wall as we can see in the Figures 3 and4d. Hence considering the expression of N u (Eq. ( 41) ) it is logical that N u increases with α * T . Similar results have been reported in [START_REF] Hong | Analysis of thermal dispersion effect on vertical-plate natural convection in porous media[END_REF] where authors showed that when the transverse dispersion effect dominates, the heat transfer is greatly increased. This is also in agreement with the results reported by [START_REF] Sheremet | Effect of thermal dispersion on transient natural convection in a wavy-walled porous cavity filled with a nanofluid: Tiwari and das nanofluid model[END_REF] for natural convection in a porous cavity filled with a nanofluid. Fig. 6b shows that N u decreases with the increase of the permeability anisotropy ratio (r k ). This is consistent with the results obtained in [START_REF] Bennacer | Double diffusive convection in a vertical enclosure filled with anisotropic porous media[END_REF] and [START_REF] Ni | Natural convection in a vertical enclosure filled with anisotropic porous media[END_REF] for natural convection (without thermal dispersion) and for equivalent ranges of parameters. This behavior can be explained by the fact that at a constant value of the Rayleigh number (i.e. k y is constant), the increase of r k can be interpreted as a decrease of the permeability in the horizontal direction k x . This entails a weaker convective flow with more expanded thermal boundary layers ( in the lower part of the hot wall) and by consequence smaller N u. The marginal effect of α * L (Fig. 6c) indicates the existence of two regimes for the evolution of the Nusselt number. Thus, we have a decreasing N u for α * L < 0.25 and an inverse behavior for α * L > 0.25. Indeed, when α * L is increased, the mixing zone between the hot and cold fluids expands in the direction of the flow. This will push the highest and lowest isotherms towards the vertical walls and increase by consequence the thermal gradient in the vertical boundary layers. On the other hand, this redistribution of the thermal gradient leads to an attenuation of the rotating flow within the cavity. Thus, referring to the expression of the Nusselt number (Eq. ( 41)), we can deduce that, for small values of α * L (< 0.25), N u decreases as the velocity variation is predominating. For large values of α * L (> 0.25), N u increases because the effect of the thermal gradient becomes significant and more important than the velocity. Fig. 7 shows bar-plots of the first order and total Sobol' indices of the maximum velocity u * max . Results indicate that the variability of u * max is mainly controlled by Ra and r k . Interactions between Ra and r k are also observed. They explain 14.5% of the total variance of u * max . The total effect of α * T accounts for approximately 1.0%. In Fig. 8, we display the marginal effect of the uncertain parameters on u * max . One can observe that u * max increases with the increase of Rayleigh number Ra, indicating that the buoyancy-induced flow along the horizontal surfaces becomes much stronger as Ra is increased. On the contrary, we note that u * max decreases with the increase of r k , as the latter corresponds to the decrease of the permeability in the horizontal direction. We can also note that small variations of u * L is associated with a decrease of u * max . The effect of α * L on the velocity can be understood with the help of the stream function form of the flow equation. This form can be obtained by applying the curl operator on the Darcy's law. It is a Poisson equation with the horizontal component of the temperature gradient as source term. This equation is subject to zero stream function as boundary conditions. The corresponding solution is concentric streamlines. The shape of these streamlines (center, orientation, spacing and density) depends on the source function. For the problem of natural convection in square cavity, the maximum horizontal temperature gradient is located at the right bottom and left top corners. The resulting streamlines have a concentric ellipsoidal shape with focal axis oriented in the direction of the line connecting the maximum gradient points (the cavity first bisector). The increase of α * L leads to the enhancement of the heat mixing by longitudinal dispersion in the zones where the velocity is parallel to the temperature gradient (outside the boundary layers of the hot and cold walls).

By consequence, the horizontal temperature gradient decreases in this zone and increases outsides. This implies that, on the one hand, the reorientation of the focal axis of the ellipsoidal shaped streamlines in the horizontal direction and, on the other hand, an attenuation of the stream function maximum value. Hence, the rotating flow decelerates, u * max decreases and the point of maximum velocity moves toward the center of the cavity surfaces. Reverse behavior can be observed when α * T is increased. In such a case the transverse heat mixing is enhanced in the zone where the velocity is orthogonal to the temperature gradient (within the boundary layers). Transverse heat flux is horizontal in this case.

Hence, the horizontal temperature gradient decreases within the boundary layers and increases outside.

The direction of the focal axis moves toward the first bisector and the maximum value of the stream function increases. The streamlines becomes more spaced at the vertical walls and more closed at the top surfaces. The value of u * max increases and its location at the top (resp. bottom) surface moves toward the cold (resp. hot) wall. Fig. 9 shows the bar-plots of the first order and total Sobol indices of the maximum velocity v * max . Results shows that the variability of v * max is mainly due to the principal effects of Ra. Unlike u * max , v * max is only slightly sensitivity to r k .

The differences between the first-order and total indices are negligible, indicating insignificant interactions between the parameters. Fig. 10 represents the marginal effect of the different parameters on the v * max . As expected, a different magnitude of variations is obtained indicating the level of influence of the parameters. The largest variation is obtained with the Rayleigh number Ra. Fig. 10a shows that v * max increases with the increase of Ra as this latter intensify the rotating flow within the cavity. The marginal effect of v * max to the permeability ratio r k is slightly flatter (see Fig. 10b) and confirm the weak sensitivity of v * max to r k . The obtained negative slope is the consequence of the horizontal velocity reduction caused by the decrease of k x . This finding is consistent with the results obtained in [START_REF] Bennacer | Double diffusive convection in a vertical enclosure filled with anisotropic porous media[END_REF]. As expected, the marginal effects of v * max to α * L and α * T are nearly flat (Fig. 10(c-d

Uncertainty quantification

The constructed PCE can be employed as a surrogate model of the target output. Statistical analysis is then affordable by performing a large Monte Carlo simulations on the PCE approximations at a very low additional computational cost upon sampling the random input parameter space. We depict in Fig. 11 the probability density function (PDF) of the scalar QoIs. The PDFs computed by relying on the PCE using 10 5 Monte Carlo simulations are compared with the PDFs computed with the 1, 000 Monte Carlo (MC) simulations of the physical model. Note that we limit our comparison to 1, 000

MCs of the complete model because of the computational burden. A number of conclusions can be drawn with respect to Fig. 11. The marginal PDFs resulting from the PCEs compare very well with the PDF obtained by relying on numerical MC simulations at a total cost of 150 simulations though.

Positively skewed distributions, with longer tails toward larger values are observed for both output N u and u * max . The long tail of the PDF of the N u is associated with setting characterized by low values of α * T , whereas the long tail of u * max is associated with setting characterized by low values of r k . A flat distribution with short tails is obtained for v * max . This is due to the fact that v * max is mainly sensitive to the Rayleigh number Ra and that they are linearly related.

In the following, we investigate the level of correlation between QoI's pairs by determining the correlation coefficient, defined as the covariance of the two output variables divided by the product of their standard deviation. Table 3(a) lists the correlation coefficient evaluated using 1, 000 MC simulations. To further elaborate on the accuracy of the PCE, the correlation coefficient results obtained by relying on PCE are also shown (see Table 3(b)). Again, the agreement between the full model and the PCE is quite remarkable. A strong positive correlation between u * max and v * max is observed. This is related to the fact that u * max and v * max are both affected by the Rayleigh number Ra and they both increase with its increase.

Effect of heterogeneity

In this section, an heterogeneous porous media will be considered. The heterogeneity is assumed to follow the exponential model given in Eq. ( 9) and ( 10). The effect of heterogeneity is expressed in terms of the rate of change of the permeability σ * . Consequently, five independent input random variables X = {Ra, r k , α * L , α * T , σ * } are now considered uncertain with uniform marginal distributions.

The spatial discretization required to reach the converged numerical solution is highly dependent on the degree of heterogeneity. Indeed, an increase in the heterogeneity degree results in an increase of the local Rayleigh number, leading to a locally steeper and rougher solution than for the homogeneous case Shao et al. (2016). Consequently, finer meshes are required to obtain a mesh independent solution. As for the homogeneous case, the most challenging configurations of parameters are thoroughly tested.

A special irregular mesh is used to obtain the converged finite element solution. This mesh involves local refinement on the high-permeability zones where the buoyancy effects are more significant. A non uniform grid of 64, 000 nodes is used. All simulations were performed for 8[t] in order to be sure that the steady state solution is reached. These discretization parameters are kept fixed in subsequent simulations.

In view of computing the PCE expansion of the model outputs, an experimental design drawn with QMC of size N = 150 is considered. As in the previous case, the candidate basis is determined using a standard truncation scheme with q = 1 for all the outputs. The corresponding results (e.g. polynomial degree giving the best accuracy, relative LOO error and number of retained polynomials) of the PCE are given in The corresponding depictions of the temperature variance is shown in Fig. 12b. First, we observe is that the spatial distribution of the variance is not anymore symmetric as in the homogeneous case.

The heterogeneity results in a different distribution of the variance while maintaining the same level of variability. The smallest variation zone is located at the boundary layer of the vertical walls and in the slow-motion region, as for the homogeneous case. It should be noted here that, due to the effect of increased permeability, the slow-motion region expands horizontally and moves up toward the right corner Fahs et al. (2015a). The zone of low temperature variance exhibits similar behavior. The largest variance zone moves to the low permeable layers and it is shifted toward the hot wall. The high permeability at the top layers leads to a reduction of the variance of the temperature.

The spatial maps of the total Sobol' indices are depicted in Fig. 13. It demonstrates that, due to the heterogeneity, the symmetry of the Sobol' indices spatial distribution around the center of the cavity is completely destroyed. Fig. 13 indicates that α * L , α * T and σ * have significant influence on the temperature distribution.

A closer look to the temperature variance confirms that, as for the homogeneous case, α * L is the most sensitive parameter since its zone of influence intersects well with the zone of maximum temperature variance. Comparing to the homogeneous case, the zone of influence of α * L expands in the zones of low in the lower part of the cavity. The reverse is true for the zone of influence of α * T . This zone expands vertically near the hot wall and contracts near the cold wall. This behavior is also attributable to the shifting of the slow-motion region toward the top right corner due to heterogeneity. The sensitivity of the temperature distribution to the rate of change of heterogeneity σ * is mainly important around the zone of slow rotating motion as it can be seen in Fig. 13e. The zone of influence of Ra expands around the right bottom corner (Fig. 13a). This behavior is related to the expansion of the zone in which the velocity is relatively weak as a consequence of the low permeability near the cavity bottom surface. Conversely, around the top left corner the high permeability induces faster convective flow associated with higher dispersion tensor. Hence, mixing by thermal dispersion dominates and temperature distribution becomes less sensitive to Ra. Fig. 13b shows that in the case of heterogeneous porous media, r k becomes more influential on the temperature distribution than for the homogeneous case. Its zone of influence expands around the top left corner.

• GSA of the scalar QoIs • The heterogeneity of the porous media does not affect the rank of the parameters regarding their influences on the model outputs. Ra and α * T remain the most influential parameters for N u, Ra and r k for u * max , and Ra for v * max .

• The uncertainty associated with the rate of change of the heterogeneity σ * has no effect on N u. This is coherent with the results obtained for the temperature distribution that shows that the effect of σ * is located relatively far from the hot wall in the slow-motion region.

• The influence of σ * is more important on u * max and v * max , which is reasonable because the velocity is directly related to the permeability level of heterogeneity.

• One can also observe that the heterogeneity renders N u and v * max more sensitive to r k . In fact, in vertically stratified heterogeneous domains (as it is the case here), the anisotropy can be at the origin of a vertical flow that can transmit the fluid from a certain layer to lower or upper ones with different permeability. Hence any change of r k may have a strong effect on the flow as different layers of heterogeneity can be involved.

The investigation of the marginal effects of the parameters on the model outputs showed that the conclusions drawn for the homogeneous case still hold for the heterogeneous porous media. This is why we hereby discuss only the marginal effect of σ * . The results show a slight variation of N u against σ * . Indeed, the increase of sigma leads to the attenuation of the local velocity at the lower part of the hot wall and an intensification in the upper part. The attenuation of the horizontal velocity around the bottom surface is associated to a reduction of the thermal dispersion and by consequence a decrease of the temperature gradient. Reverse process occurs at the upper surface and leads to the increase of the thermal gradient. As a consequence, the local Nusselt number decreases in the lower part of the

Summary and conclusions

The proposed study handled the topic of natural convection problem in heterogeneous porous media including velocity-dependent dispersion. The considered benchmark is the popular square cavity filled with a saturated porous medium and subject to differentially heated vertical walls and adiabatic horizontal surfaces. The simplicity of the geometry and the boundary conditions renders this problem especially suitable for testing numerical models but also useful to provide physical insights into the involving processes. It introduces however several uncertain parameters, namely: the average Rayleigh number (Ra), the permeability anisotropy ratio (r k ), the non-dimensional dispersion coefficients (α * L and α * T ).

The imperfect knowledge of the system parameters and their variability significantly affect the flow and heat transfer patterns. It is thus of utmost importance to properly account for the aforementioned uncertainties within the frame of uncertainty and sensitivity Analysis. In this work, we analyze the impact of the uncertain parameters on the output quantities of interest (QoIs) allowing assessment of flow and heat transfer. To describe the flow process, we use the maximum dimensionless velocity components (u * max and v * max ). For the heat transfer process, the assessment is based on the spatial distribution of the dimensionless temperature (T * ) and the average Nusselt number N u on the hot wall. The effect of heterogeneity was also investigated by considering a stratified heterogeneous porous media with an exponential distribution of the permeability as a function of depth. The rate of change of the permeability σ * is then considered uncertain.

Herein, we performed a comprehensive global sensitivity analysis and uncertainty quantification through performing the probability distributions by means of surrogate modeling. Sparse PCE are used for this purpose. Our results lead to the following major conclusions.

• Sparse PCE proved particularly efficient in providing reliable results at a considerably low computational costs. All results derived for the homogeneous (resp. heterogeneous) case were obtained at the cost of only 150 simulations of the computational model. Note that those runs could have been carried out in parallel on any distributed computing architecture. Due to the inexpensiveto-evaluate formulation of the sparse PCE meta-model, the probability density functions (PDF's) of the QoIs have been accurately estimated. An excellent agreement between sparse PCE and MC results is obtained.

• The Sobol' indices of the temperature distribution allow specifying the spatial zone of influence of each parameter. Results showed that the variability of the temperature distribution is largely influenced by the effect of α * L and α * T . Nevertheless, the effect of α * L on the temperature distribution is more pronounced than that for α * T as its zone of influence is located in the region where the variance is maximum.

• The variability of N u is mainly due to main effects of Ra and α * T . Indeed, the Rayleigh number dramatically influence the flow profile and heat transfer within the cavity, as well as the thermal boundary layer thickness. The variability of u * max is mainly due to main effects of the ratio of anisotropy r k , and Ra while the variability of v * max is mainly controlled by Ra. • The effect of the heterogeneity results in a different distribution of the variance of the temperature while maintaining the same level of variability. The zone of largest temperature variance becomes located in the low permeable layer near the bottom surface of the cavity. In this case, the isotherms are more affected by the circulating flow than in the homogeneous case, especially in the high permeable zones, associated to a more intense convective flow. Results shows that the average Nusselt number is not sensitive to the heterogeneity rate while an increase of σ * is associated to an increase of the maximum velocity u * max and v * max . • Marginal effects of each parameters are also readily obtained from PCE. As opposed to classical "one-at-a-time" sensitivity analyses, where all parameters but one are frozen so as to study the effect of the remaining one, the univariate effect curves account for the uncertainties in the other parameters. Quantitative conclusions have been drawn, which confirm qualitative interpretations.

This study represents a prototype to point out the benefit of GSA and UQ to understand how the complex system of natural convection in porous media behaves. This kind of study shall be useful for safe design and risk assessment of systems involving natural convection in porous enclosure.

Fig. 1 :

 1 Fig. 1: Schematic diagram of the heterogeneous porous-cavity problem

Fig. 2 :

 2 Fig. 2: Homogeneous case: Comparison between PCE and the true model on 1, 000 validation runs for the model outputs.

  Fig. 3: Homogeneous case -spatial distribution of the temperature statistical moments.

Fig. 4 :

 4 Fig. 4: Homogeneous case -Spatial distribution of the total Sobol' indices for the temperature

Fig. 5 :

 5 Fig. 5: Homogeneous case -Total Sobol' indices for the average of the average of Nusselt number N u;

Fig. 6 :

 6 Fig. 6: Homogeneous case -Univariate effects of the input parameters on the average Nusselt number N u

Fig. 7 :

 7 Fig. 7: Homogeneous case -Total Sobol' indices for u * max ;

Fig. 8 :Fig. 9 :Fig. 10 :

 8910 Fig. 8: Homogeneous case -Univariate effects of the input parameters on u * max

Fig. 11 :

 11 Fig. 11: Homogeneous case -Marginal pdf obtained by relying on PCE (red lines) and 1, 000 MC Simulations (blue lines).

Fig. 12a illustrates

  Fig.12aillustrates the spatial distribution of the mean temperature by relying on the PCE. It shows that the isotherms are more affected by the circulating flow than that of the homogeneous case; especially in the high permeable zones (near the top surface of the cavity). Indeed, the local Rayleigh number exceeds the average value Ra, leading eventually to a more intense convective flow.

Fig. 12 :

 12 Fig. 12: Heterogeneous case -spatial distribution of the temperature statistical moments.

(

  Fig. 13: Heterogeneous case -Spatial distribution of the total Sobol' indices for the temperature

Fig. 14 Fig. 14 :

 1414 Fig. 14 shows bar-plots of the first order and total Sobol' indices of the three model outputs ( N u, u * max and v * max ).

  vertical wall and increases at the upper part. Upper and lower local Nusselt numbers tend to balance out each other and lead to small variation of the N u. Marginal effects of σ * on u * max and v * max are given in Fig. 15. This figure confirms the high sensitivity of u * max and v * max to sigma. It indicates also an increasing variation of u * max and v * max against σ * as a result of the permeability increasing at the top surface of the cavity.

Fig. 15 :

 15 Fig. 15: Heterogeneous case -Univariate effects of σ * on u * max and v * max .

  In this particular context, the Sobol sensitivity indices Sobol' (1993);[START_REF] Homma | Importance measures in global sensitivity analysis of non linear models[END_REF]; Sobol' (2001) are widely used as sensitivity metrics, because they do not rely on any assumption regarding the linearity or monotonous behavior of the physical model. Various techniques have been proposed in the literature for computing the Sobol indices, see e.g.[START_REF] Archer | Sensitivity measures, ANOVA-like techniques and the use of bootstrap[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF];[START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF];

Table 1 :

 1 Results of the utilized PCE

		N u	u * max	v * max
	p opt	5	6	4
	err LOO Size of the Sparse Basis	8.6 × 10 -4 5.81 × 10 -2 3.93 × 10 -4 63 47 41
	different behaviour of the parameters is observed. The first and total Sobol' indices are computed, as

well as second-order one, based on the obtained PCEs of the various QoI's. Referring to the results in Table

Table 2 :

 2 Homogeneous case -Sobol' indices for N u

	Ra	r k	α * L	α * T
	S i 0.739 0.0547 0.0252 0.1093
	S T 0.803 0.077 0.0235 0.158

Table 3 :

 3 Correlation coefficient of the QoIs evaluated with 1, 000 MC simulations (Black fonts) and PCE.

  Table 4 for the three scalar output N u, u * max and v * max . An accurate PCE is obtained for both N u and v * max , where LOO error is about 1%. A less accurate PCE is obtained for u * max , where LOO error is larger than 0.1.

Table 4 :

 4 Results of the utilized PCE

		N u	u * max	v * max
	p opt	6	5	6
	err LOO size of the sparse basis	1.1 × 10 -2 2.88 × 10 -1 9.68 × 10 -3 46 26 63
	• GSA of the temperature field


Guillaume Devailly

Katia Fêve

Safia Saci

Julien Sarry

Sophie Valière

Olivier Bouchez

Laure Ravon

Yvon Billon

Martin Beaumont

Hélène Gilbert

Transcriptomic and DNA methylation response to feed intake in the duodenum in high-and low-feed efficiency pig lines

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Guillaume Devailly 1 , Katia Fêve 1 , Safia Saci 1 , Julien Sarry 1 , Sophie Valière 2 , Olivier Bouchez 2 , Laure Ravon 3 , Yvon Billon 3 , Martin Beaumont 1 , and Hélène Gilbert 1 1: GenPhySE, Université de Toulouse, INRAE, ENVT, 31326, Castanet Tolosan, France 2: INRAE, US 1426, GeT-PlaGe, Genotoul, Castanet-Tolosan, France 3: Pig phenotyping and Innovative breeding facility, GenESI, UE1372, INRAE, F-17700 Surgères, France

Key Results

• 3 times as many genes differentially expressed due to feed intake than between the two pig lines. • The duodenum transcriptomic response to feed intake is greater in the more efficient pig line. • Duodenum DNA methylation profiles are barely affected by short-term feed intake. • Duodenum DNA methylation profiles are distinct between the two pig lines.

Futur directions

• Integration of transcriptomic and DNA methylation profiles to investigate if differences in DNA methylation might explain part of the different transcriptomic response between pig lines. • Transcriptomic and DNA methylation profiling of the stomach mucosa.

• Cleaner identification of differentially methylated regions by focusing on unmethylated regions.

This project is funded by the department of Animal Genetics, INRAE. RNA-seq and MeDP-seq data will be publicly available soon at the European Nucleotide Archive (PRJEB46060)

Occitanie -Toulouse, France

Transcriptomic and DNA methylation response to feed intake in the duodenum in high-and low-feed efficiency pig lines Context Feed efficiency is a complex trait of great interest in animal breeding, notably to reduce environmental impact of livestock herding and to reduce the cost of animal feed for the farmers. A successful divergent selection on feed efficiency was carried out in pigs for more than 10 generations at an INRAE experimental unit, establishing lines of related good and poor feed efficiency large white pigs. The objective of the present study was to identify the molecular mechanisms underlying the divergence in feed efficiency between the two pig lines at the intestinal level.

Method

Duodenum mucosal samples were collected in post-weaning pigs from the two divergent lines, in two groups: one after a 12 hours period of feed restriction, the other with a 10 hours period of feed restriction followed by a 2 hours period of ad libitum feed access (n=6 per group and per line). Transcritome and methylome of duodenal mucosa were analyzed by poly-A RNA sequencing (RNA-seq) and by methylated DNA precipitation followed by sequencing (MeDP-seq), respectively.

> Sample production > DNA methylation analysis

By methylated-DNA precipitation followed by high throughput sequencing (MeDP-seq). Differentially methylated regions (DMR) identified by four different methods:

  11 th generation: 60 days-old G11+ (HRFI, poor feed efficiency), and G11-(LRFI, good feed efficiency).Duodenum mucosa dissectionBalanced experimental groups:• 24 pigs in total • 12 G11+ and 12 G11-•12 fasted and 12 fed • 13 castrated males and 11 females • 4 pigs each from 3 litters per line Focus on the Incretin pathway /rnaseq Mapping and quantification: Salmon DGE: limma-voom Expression ~ Condition + Line + Sex This experimentation was authorized by the French Minister of Higher Education, Research and Innovation (APAFIS#21107-2018120415595562 v10), after examinations by the animal experimentation ethic committee 084.No significant category was detected in genes over expressed in G11-.