Hossein Davarzani
email: h.davarzani@brgm.fr

Manuel Marcoux

& Michel Quintard

A LOCAL THERMAL NON-EQUILIBRIUM MODEL FOR COUPLED HEAT AND MASS TRANSFER WITH DISPERSION AND THERMAL DIFFUSION IN POROUS MEDIA

A theoretical model is proposed to predict Darcy-scale mass transport in porous media coupled with non-equilibrium heat transfer and taking into account the thermal diffusion process. A volume-averaging technique was used with approximations leading to a two-equation or two-temperature model for the macro-scale energy balance equations. Because of the Soret effect, the concentration deviation with respect to the averaged value is a function of concentration and temperature gradients in the fluid phase, the temperature gradient in the solid phase, and the difference between the two averaged temperatures. The mapping between deviations and averages involves four closure problems for the mass transport equations: problems which were solved numerically over a two-dimensional periodic-unit cell for evaluation purposes. The results show that the effective coefficients depend strongly on the thermo-physical properties of the medium and the Péclet number. In particular, the effective-Soret coefficient in porous media changes with the Péclet number and the phases' thermal conductivity ratios.

In many natural and industrial processes we must take into account the coupling between mass transport and temperature. When this occurs, an anti-diffusive mass flux can be generated in a fluid mixture under the existing temperature gradient, which tends to separate the species. This cross-effect phenomenon is known as the Soret effect or thermal diffusion. The term thermal diffusion, first used for gas mixtures, later came to refer to the phenomenon for all phase states.

Thermal diffusion can occur in many processes involving coupled heat and mass transfer in porous media, for instance, in the chemical industry [START_REF] Köhler | The Soret Effect in Liquid Mixtures -A Review[END_REF][START_REF] Reddy | Thermal Diffusion and Diffusion Thermo Impact on Chemical Reacted MHD Free Convection from an Impulsively Started Infinite Vertical Plate Embedded in a Porous Medium Using FEM[END_REF], contaminant fate and transport and environmental processes [START_REF] Beyer | Simulation of Temperature Effects on Groundwater Flow, Contaminant Dissolution, Transport and Biodegradation Due to Shallow Geothermal Use[END_REF][START_REF] Mon | Temperature Effects on Solute Diffusion and Adsorption in Differently Compacted Kaolin Clay[END_REF][START_REF] Perrier | Seasonal Thermal Signatures of Heat Transfer by Water Exchange in an Underground Vault[END_REF][START_REF] Ratchagar | Contaminant Transport Modeling Due to Thermal Diffusion Effects with the Effect of Biodegradation[END_REF], drying or vitrification processes [START_REF] Dragun | Geochemistry and Soil Chemistry Reactions Occurring during In Situ Vitrification[END_REF], and reservoir engineering [START_REF] Firoozabadi | Thermodynamics of Hydrocarbon Reservoirs[END_REF][START_REF] Galliéro | On Thermal Diffusion in Binary and Ternary Lennard-Jones Mixtures by Non-Equilibrium Molecular Dynamics[END_REF][START_REF] Vaerenbergh | Multicomponent Transport Studies of Crude Oils and Asphaltenes in DSC Program[END_REF].

In a previous study (Davarzani et al., 2010a), a theoretical model was developed using the volume-averaging technique to predict coupled heat and mass transfer equations in porous media at the Darcy scale, in the case of local-equilibrium heat transfer (i.e., with a single temperature for all phases). The results revealed that in the macro-scale model the Soret coefficient, defined as the ratio of the thermal diffusion coefficient to the binary-diffusion coefficient, can be different from the free-liquid coefficients when convection is dominant. However, the Soret coefficient was found to be the same in a free fluid and in a porous medium for a diffusiondominated regime (Davarzani et al., 2010b;[START_REF] Platten | The Soret Coefficient in Porous Media[END_REF]. In all cases, the Soret coefficient in porous media is independent of the solid phase particle-particle contact or connectivity [START_REF] Davarzani | Effect of Solid Thermal Conductivity and Particle-Particle Contact on Effective Thermodiffusion Coefficient in Porous Media[END_REF].

In a convection-dominant regime, the thermal conductivity ratio between fluid and solid plays an important role in determining the effective thermal diffusion coefficients (Davarzani et al., 2010a). The proposed theoretical model has been validated with experimental results for pure diffusion cases. A modified two-bulb method has been used to evaluate the free flow and effective diffusion as well as thermal diffusion coefficients in granular porous media (Davarzani et al., 2010b). This validated theoretical model can capture the coupled heat-and mass-transfer behaviors when the medium's thermo-physical contrasts are moderate, a situation that favors local equilibrium. However, when the thermo-physical contrast becomes high, the situation is much more complex and may lead to non-local thermal-equilibrium conditions, which is the focus of this paper. Davarzani et al. (2010) compared averaged temperatures predicted from an upscaled macroscopic model and the results coming from the solution of microscopic pore-scale equations over a periodic array of disks. The results did indeed show that the theoretical local-equilibrium model can predict transport phenomena efficiently for moderate thermal-property contrasts. In this case, the results are less dependent on the boundary and initial conditions. For higher thermo-physical contrasts, they observed a significant difference between predicted temperature and the pore-scale averaged temperatures (solid and liquid phases) in the transient period. By contrast, at steady-state and for the considered initial boundary-value problem, the theoretical local-equilibrium model can again predict how the system will behave correctly (Davarzani et al., 2010a). There is therefore an essential need to develop a local non-equilibrium model if the interest is in transient behaviors with higher phase-contrast properties. The main objective of this work was to develop such a complex and advanced model.

It is important to recognize that non-equilibrium conditions cannot be described in general by a single macro-scale model. Therefore, alternative models have been developed to take into account the local non-equilibrium problem. There are at least two different approaches to treat this problem: the two-equation model (also called dual-porosity or two-region model) and the hybrid (or mixed) model. These models are well known in petroleum engineering for the description of the fractured mediums, in chemical engineering, and hydrogeology (mobilemobile and mobile-immobile systems). The two-equation model consists of two coupled Darcyscale balance equations. These two equations are connected through an exchange term. The domains of validity of different models, including two-equation, one-equation local equilibrium, and one-equation time-asymptotic non-equilibrium models, which depend on the Péclet number and a characteristic time, have been explored by [START_REF] Davit | Modeling Non-Equilibrium Mass Transport in Biologically Reactive Porous Media[END_REF]Davit and Quintard, 2015a).

While the one-equation local equilibrium model and the one-equation non-equilibrium asymptotic model are of great interest, it was shown in [START_REF] Davit | Modeling Non-Equilibrium Mass Transport in Biologically Reactive Porous Media[END_REF] that they have limitations. In particular, it is impossible with these models to assess whether non-equilibrium conditions are or have been reached. It is clear that the local equilibrium model fails to represent actual fields in non-equilibrium conditions. Asymptotic conditions may have been reached, thus providing some justification for the use of the asymptotic model. However, these conditions cannot be verified with this model. In addition, if the quasi-steady solution reverts to local equilibrium, the asymptotic model yields too large an effective-thermal tensor and thus fails to capture actual fields. An example of such situations is given in Davarzani et al. ( 2010) (Davarzani et al., 2010a).

Contrary to the traditional two-equation model where both equations for fluid and solid are written at the macro-scale, i.e., they come from the volume averaging, the mixed models couple the macroscopic description of one continuous phase for fluids with a pore-scale equation for the solid particles [START_REF] Arbogast | Analysis of the Simulation of Single Phase Flow through a Naturally Fractured Reservoir[END_REF][START_REF] De Swaan | Analytic Solutions for Determining Naturally Fractured Reservoir Properties by Well Testing[END_REF][START_REF] Kazemi | Pressure Transient Analysis of Naturally Fractured Reservoirs with Uniform Fracture Distribution[END_REF]. It was found that the mixed model matches better the direct numerical simulation results compared to the two-equation model when the thermal conductivity of the solid phase is orders of magnitude lower than the thermal conductivity of the flowing fluid (there is also a condition often overlooked in the literature that volume fractions are not close to unity or zero. For instance, in the framework of asymptotic homogenization, it must be on the order of ε 2 (ε being in the context of the lengthscale ratio between the pore-scale and Darcy-scale). Recently, [START_REF] Xin | Numerical Investigation of Local Thermal Non-Equilibrium Effects in Coal Porous Media with Cryogenic Nitrogen Injection[END_REF] studied a mixed model for high permeability granular porous media with very low solid-phase thermal conductivity. As expected, they showed that the equilibrium time of the fluid-solid center is longer for the low thermal conductivity of the solid matrix. This equilibrium time is affected significantly by solid size and flow rate [START_REF] Xin | Numerical Investigation of Local Thermal Non-Equilibrium Effects in Coal Porous Media with Cryogenic Nitrogen Injection[END_REF][START_REF] Xin | Numerical Investigation of Local Thermal Non-Equilibrium Effects in Coal Porous Media with Cryogenic Nitrogen Injection[END_REF]. Mixed model validity is limited to high thermos-physical contrasts of phases (e.g., very low solid/fluid thermal conductivity ratios) encountered in mobile-immobile and dispersed solid-phase systems [START_REF] Debenest | Transport in Highly Heterogeneous Porous Media: From Direct Simulation to Macro-Scale Two-Equation Models or Mixed Models[END_REF][START_REF] Xin | Numerical Investigation of Local Thermal Non-Equilibrium Effects in Coal Porous Media with Cryogenic Nitrogen Injection[END_REF]. Moreover, in the case where the local nonequilibrium situations come from the fluid phase (e.g., heat exchangers, impact of heat sources), the mixed model does not provide a suitable solution. Also, the computational time for a mixed model is relatively high because it involves pore-scale grids and can become close to the calculation time of the direct numerical simulations if more precision is needed [START_REF] Debenest | Transport in Highly Heterogeneous Porous Media: From Direct Simulation to Macro-Scale Two-Equation Models or Mixed Models[END_REF] and if the population of grains inside the averaging volume cannot be replaced by the dynamics of a single grain or another type of simplified volume. Therefore, special attention should be paid to select a suitable local non-equilibrium model according to the problem and its context.

While the development of the two-equation model is more complicated, in many applications it gives better and more robust results. The lower required computational time is also a good argument to use this model in many large-scale applications such as soil contamination and petroleum engineering. It also must be emphasized that transient closures, involving time convolution with effective properties coming from the unsteady-state version of the closure problems, provide an accurate transient behavior that does not suffer from the limitations of mixed models, while the implementation of time-convolution in the macro-scale numerical model is not straightforward (see [START_REF] Chastanet | Mass Transfer Process in a Two-Region Medium[END_REF]Davit and Quintard, 2015b)).

Finally, transient closures can be used as a basis to generate multi-rate or N-equations models as discussed in [START_REF] Davit | A Domain Decomposition Approach to Finite-Epsilon Homogenization of Scalar Transport in Porous Media[END_REF]Davit and Quintard, 2015b;[START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF], which provide an alternative to transient closure without the numerical implementation problems. Therefore, the most versatile model and easily implemented, valid but with a variable degree of accuracy for all range of parameters, is the two-equation model with steady-state closure. This must be the first choice when working on a system expecting a non-equilibrium behavior. Hence our motivation to develop such a model in the presence of thermal diffusion.

So far, no theoretical predictive model has been proposed to simulate macro-scale mass transfer with thermal diffusion in porous media, in the case of local thermal non-equilibrium. We propose such a model here by extending previous results obtained using a volume-averaging technique in the case of a local non-equilibrium heat-transfer situation with no Soret effect. First, we present the pore-scale equations with an introduction to the upscaling technique we used.

Then we recall the results (macro-scale equations, closure problems, and effective thermal conductivities) for the macro-scale energy equation in the case of local equilibrium and local non-equilibrium. Using these results, we develop the macro-scale equation for the heat and mass transport equations. Finally, we calculate the effective properties for simple representative-unit cells in order to investigate the impact of the pore-scale physical parameters on the effective properties.

GOVERNING MICROSCOPIC EQUATION

The problem under consideration corresponds to non-isothermal mass transport in a porous medium as represented schematically in Fig. 1. A binary-fluid mixture flows through a saturated porous medium subjected to a temperature gradient. For simplicity, we follow the same notation as in [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF] and (Davarzani et al., 2010a), where the fluid phase is represented by the β-phase and the porous medium solid phase is identified as the σ-phase.

The pore-scale energy balances for the fluid (β-phase) and solid (σ-phase) with their associated boundary conditions are introduced as

�𝜌𝑐 𝑝 � 𝛽 𝜕𝑇 𝛽 𝜕𝜕 + �𝜌𝑐 𝑝 � 𝛽 𝑣 𝛽 . 𝛻𝑇 𝛽 = 𝛻. �𝑘 𝛽 𝛻𝑇 𝛽 �, in the β-phase (1)
BC1: In this study, thermal diffusion (or the Soret effect), i.e., mass transport caused by the thermal gradient, is taken into account. Therefore, the pore-scale mass balance equation and the interface-boundary condition can be presented as follows [START_REF] Nield | Convection in Porous Media[END_REF] ( ) ( )

𝑇 𝛽 = 𝑇 𝜎 , at βσ A (2) BC2: ( ) ( ) σ σ βσ β β βσ T k T k ∇ = ∇ . . n n , at βσ A (3) ( ) ( ) σ σ σ σ ρ T k t T c p ∇ ∇ = ∂ ∂ . , in the σ-phase (
β β β β β β β T D c D c t c T ∇ + ∇ ∇ = ∇ + ∂ ∂ . . v
, in the β-phase (5) BC1:

( ) 0 . = ∇ + ∇ β β β β βσ T D c D T n , at βσ A (6)
where β c is the species mass concentration in the β-phase, β D is the binary molecular diffusion, and β T D is the thermal diffusion coefficient in the free fluid. Here, we have considered that the fluid-solid interface is inert and impervious. Therefore, we have neglected any solute accumulation and exchange phenomena (e.g., adsorption, dissolution, reaction) at the β-σ interface.

We have not taken into account the Dufour effect in the energy equation here. This is generally an acceptable approximation. It has the fundamental advantage of decoupling the heatand mass-transport problems, provided that fluid properties do not also depend on the concentration. Consequently, we treated the upscaling problem in a sequential manner:

temperature first (This problem has received a considerable attention in the literature [START_REF] Carbonell | Heat and Mass Transfer in Porous Media, Fundamentals of Transport Phenomena in Porous Media[END_REF][START_REF] Levec | Longitudinal and Lateral Thermal Dispersion in Packed Beds. Part I: Theory[END_REF][START_REF] Quintard | Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF][START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Vortmeyer | Equivalence of One-and Two-Phase Models for Heat Transfer Processes in Packed Beds: One Dimensional Theory[END_REF][START_REF] Zanotti | Development of Transport Equations for Multiphase System-I: General Development for Two Phase System[END_REF]) and then the mass-transfer problem, in which heat-transfer closure was already known.

In this study, we have assumed that density and viscosity are constant and independent of temperature and concentration, such that the momentum equations are not coupled with the energy and mass-conservation equations. In fact, it is sufficient that this approximation holds at the REV-scale. In this case, variations of the fluid properties can be taken into account in the macro-scale equations. Therefore, the fluid velocity fields are presumed to be known beforehand when dealing with those equations.

UPSCALING TECHNIQUE

There are usually computational limitations to simulating transport phenomena at the pore scale in a large porous-medium domain, due to the required computational memory and time. In practice, a homogenous-equivalent model can simplify modeling multi-phase transport in porous media [START_REF] Brenner | Dispersion Resulting from Flow through Spatially Periodic Porous Media[END_REF][START_REF] Carbonell | Heat and Mass Transfer in Porous Media, Fundamentals of Transport Phenomena in Porous Media[END_REF][START_REF] Cushman | A Primer on Upscaling Tools for Porous Media[END_REF][START_REF] Mei | Method of Homogenization Applied to Dispersion in Porous Media[END_REF]. Several theoretical methods and techniques exist for developing such a model and moving through different scales [START_REF] Cushman | A Primer on Upscaling Tools for Porous Media[END_REF]. The volume-averaging method [START_REF] Carbonell | Heat and Mass Transfer in Porous Media, Fundamentals of Transport Phenomena in Porous Media[END_REF] and homogenization techniques [START_REF] Mei | Method of Homogenization Applied to Dispersion in Porous Media[END_REF] are the most widely used techniques. They share some common principles, in particular the idea of looking at the pore-scale field in terms of a correction or deviation to the macro-scale field, as discussed in [START_REF] Davit | Homogenization via Formal Multiscale Asymptotics and Volume Averaging, How Do the Two Techniques Compare?[END_REF]. Our paper is an extension of previous papers devoted to the development of local non-equilibrium heat-transfer models using the framework of a volume-averaging theory. We will follow these developments closely and emphasize the original points when the Soret effect plays a role. The volume-averaging technique consists of first producing the macro-scale equations from averaging the pore-scale equations; second, solving micro-scale problems over a unit cell to obtain the effective properties; and finally, solving the upscaled macro-scale model with calculated effective properties and specific initial and boundary conditions.

Our starting point in this paper, given the above arguments on the sequential coupling between the heat-and mass-transport problems, is the local non-equilibrium two-temperature model extensively studied in the literature, especially using the volume-averaging method [START_REF] Degroot | Closure of Non-Equilibrium Volume-Averaged Energy Equations in High-Conductivity Porous Media[END_REF][START_REF] Puiroux | Non-Equilibrium Theories for Macroscale Heat Transfer: Ablative Composite Layer Systems[END_REF][START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Quintard | One-and Two-Equation Models for Transient Diffusion Processes in Two-Phase Systems[END_REF]. The method and its application to the heat-transfer problem will be described without all the mathematical developments, as a preliminary for the original part of the paper: the treatment of the mass-transport problem with the Soret effect and thermal non-equilibrium.

To obtain a Darcy-scale model for the coupled heat-and mass-transfer problem, the porescale equations must be integrated over a representative elementary volume, V shown in Fig. 1.

While a more general and accurate averaging procedure may be used ( [START_REF] Davit | Technical Notes on Volume Averaging in Porous Media I: How to Choose a Spatial Averaging Operator for Periodic and Quasiperiodic Structures[END_REF]Quintard and Whitaker, 1994a;Quintard and Whitaker, 1994b;Quintard and Whitaker, 1994c;[START_REF] Quintard | Transport in Ordered and Disordered Porous Media Iv: Computer Generated Porous Media for Three-Dimensional Systems[END_REF]Quintard and Whitaker, 1994e)), the superficial and intrinsic averages associated with the phases are in the classical sense respectively defined for a field variable 𝜙 𝛽 as [START_REF] Whitaker | The Method of Volume Averaging[END_REF])

dV V ∫ = β β β ϕ ϕ V 1 (7) dV V V ∫ = β β β β β ϕ ϕ 1 (8)
where β V is the volume of the β-phase in the averaging volume. The superficial and intrinsic averages are related by

β β β β ϕ ε ϕ = (9)
where β ε is the β-phase volume fraction that will be considered as a constant in this study.

The sum of the volume fractions of the two phases satisfies

1 = + σ β ε ε (10)
Pore-scale deviations are defined as the difference between the point microscopic values and intrinsic averaged values. These spatial decompositions can be defined in the β-phase and σ -phase as

β β β β ϕ ϕ ϕ + = and σ σ σ σ ϕ ϕ ϕ + = (11)
As usual, we will make use of the following length-scale constraints.

<< << (12)

This hypothesis leads to many simplifications at different steps of the upscaling processes.

Below the macroscopic-coupled models obtained from upscaling are presented. In section 4, numerical results for the effective-diffusion and thermal-diffusion coefficients will be presented and discussed.

Macroscopic Energy-Balance Equation

We remind the reader that, with the adopted assumptions as explained in Section 2, the heattransfer problem is independent of the concentration problem [START_REF] Carbonell | Heat and Mass Transfer in Porous Media, Fundamentals of Transport Phenomena in Porous Media[END_REF][START_REF] Levec | Longitudinal and Lateral Thermal Dispersion in Packed Beds. Part I: Theory[END_REF][START_REF] Quintard | Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF][START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Vortmeyer | Equivalence of One-and Two-Phase Models for Heat Transfer Processes in Packed Beds: One Dimensional Theory[END_REF][START_REF] Zanotti | Development of Transport Equations for Multiphase System-I: General Development for Two Phase System[END_REF]. The theoretical model for temperature has received considerable attention in the literature [START_REF] Carbonell | Heat and Mass Transfer in Porous Media, Fundamentals of Transport Phenomena in Porous Media[END_REF][START_REF] Levec | Longitudinal and Lateral Thermal Dispersion in Packed Beds. Part I: Theory[END_REF][START_REF] Quintard | Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF][START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Vortmeyer | Equivalence of One-and Two-Phase Models for Heat Transfer Processes in Packed Beds: One Dimensional Theory[END_REF][START_REF] Zanotti | Development of Transport Equations for Multiphase System-I: General Development for Two Phase System[END_REF]. In this section, we review the previously described one-and two-equation models.

One-Equation Model

The Darcy-scale model for the local heat-equilibrium case is written as [START_REF] Nozad | Heat-Conduction in Multiphase Systems .1. Theory and Experiment for 2-Phase Systems[END_REF][START_REF] Quintard | Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF][START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Quintard | One-and Two-Equation Models for Transient Diffusion Processes in Two-Phase Systems[END_REF][START_REF] Whitaker | Toward a Diffusion Theory of Drying[END_REF])

( ) ( ) ( ) ( ) ( ) ( ) T T c T t c c p p p ∇ ∇ = ∇ + ∂ ∂ + . . . * k v β β β β σ σ β β ε ρ ρ ε ρ ε (13) where σ σ β β T T T = =
and * k is the effective thermal conductivity tensor defined by

β  0 r L ( ) ( ) ( ) β β β β βσ σ β σ σ β β ρ ε ε βσ T p A T c dA k k k k b v b n I k * - - + + = ∫ V (14)
where β T b is a mapping variable described below.

In this development, the spatial temperature and velocity decompositions are introduced as

β β β β T T T - = ~, in β V and σ σ σ σ T T T - = ~, in σ V (15) β β β β v v v - = ~, in β V (16)
The spatial deviations of the temperatures and average fields can be related in the upscaling process by a linear relationship as [START_REF] Carbonell | Heat and Mass Transfer in Porous Media, Fundamentals of Transport Phenomena in Porous Media[END_REF]) Finally, the closure problem with the boundary conditions at the interface of phases that must be solved numerically to find * k from Eq. ( 14) is written below (Davarzani et al., 2010a;[START_REF] Quintard | Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF][START_REF] Whitaker | The Method of Volume Averaging[END_REF] Problem I:

T T T ∇ = . ~β β b (17) T T T ∇ = . ~σ σ b (18)
( ) ( ) β β β β β β β ρ ρ T p T p k c c b v b v 2 . ∇ = + ∇ (19) BC1: σ β T T b b = , at βσ A (20) BC2: ( ) σ β βσ σ σ βσ β β βσ k k k k T T - + ∇ - = ∇ - n b n b n . . , at βσ A (21) σ σ T k b 2 0 ∇ = (22)
Periodicity:

( ) ( ) r b r b β β T i T = +  & ( ) ( ) r b r b σ σ T i T = +  , i=1,2,3 (23) 0 T T T β σ β β σ σ ε ε = + = b b b (24)

Two-Equation Model

We have not reproduced here details about the upscaling of local thermal non-equilibrium situations, i.e., situations for which the two averaged temperatures may differ significantly.

Discussion of this subject can be found in (Davarzani et al., 2010a;[START_REF] Davit | Modeling Non-Equilibrium Mass Transport in Biologically Reactive Porous Media[END_REF]Davit and Quintard, 2015a;[START_REF] Davit | Technical Notes on Volume Averaging in Porous Media I: How to Choose a Spatial Averaging Operator for Periodic and Quasiperiodic Structures[END_REF][START_REF] Davit | Correspondence between One-and Two-Equation Models for Solute Transport in Two-Region Heterogeneous Porous Media[END_REF][START_REF] Degroot | Closure of Non-Equilibrium Volume-Averaged Energy Equations in High-Conductivity Porous Media[END_REF][START_REF] Moyne | Two-Equation Model for a Diffusive Process in Porous Media Using the Volume Averaging Method with an Unsteady-State Closure[END_REF][START_REF] Quintard | Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF]. Different models may be introduced and in this paper we present the simplest and more commonly used, namely the two-temperature model [START_REF] Carbonell | Heat and Mass Transfer in Porous Media, Fundamentals of Transport Phenomena in Porous Media[END_REF][START_REF] Levec | Longitudinal and Lateral Thermal Dispersion in Packed Beds. Part I: Theory[END_REF][START_REF] Quintard | Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF][START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Vortmeyer | Equivalence of One-and Two-Phase Models for Heat Transfer Processes in Packed Beds: One Dimensional Theory[END_REF][START_REF] Zanotti | Development of Transport Equations for Multiphase System-I: General Development for Two Phase System[END_REF]. Appendix A [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF]. The resulting macroscopic equation for the β-phase can be written as

( ) ( ) ( ) ( ) σ σ β β σ σ βσ β β ββ σ σ β β β β β β β β β β β β ρ ε ρ ε T T h T T T T T c T t c T βσ ββ p p - - ∇ + ∇ ∇ = ∇ - ∇ - ∇ + ∂ ∂ . . . . . . k k u u v (27)
Here, ββ k and βσ k are the effective thermal conductivity tensors; ββ u , βσ u are the effective heat-convective vectors; T h is the heat-transfer coefficient. They can be calculated from the solution of closure problems and using the equations (A. 25), (A. 26), (A. 29), (A. 30), and (A.

24), respectively.

The intrinsic average-temperature equation for the σ-phase can be written as

( ) ( ) ( ) σ σ β β σ σ σσ β β σβ σ σ σσ β β σβ σ σ σ σ ρ ε T T h T T T T T t c T p - + ∇ + ∇ ∇ = ∇ - ∇ - ∂ ∂ . . . . . k k u u (28)
In this equation also, σβ k , σσ k , σβ u , and σσ u are effective tensors, obtained from the corresponding closure problems and expressions given in Appendix A. As mentioned in [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Quintard | One-and Two-Equation Models for Transient Diffusion Processes in Two-Phase Systems[END_REF], the effective-thermal conductivity in the case of one-equation local equilibrium model can be calculated as a limit case of the effective thermal conductivity tensors of the two-equation model [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Quintard | One-and Two-Equation Models for Transient Diffusion Processes in Two-Phase Systems[END_REF] by

σσ σβ βσ ββ k k k k k * + + + = (29)

Macroscopic Transport Equations for Concentration

The major assumptions that allow us to treat the upscaling process for the mass-transport problem in a sequential manner (by first solving the problem for β T and σ T then the problem for β c ) are:

-the Dufour effect is neglected -the pore-scale physical properties such as viscosity, density, thermal conductivity, heat capacity, and diffusion coefficients do not depend on temperature or concentration.

Consequently, β T and σ

T may be replaced in Eq. ( 5) by their averaged values and their deviations. Furthermore, those variables may be replaced by their representations, in Eqs ( 17) and ( 18), or ( 25) and ( 26). This can be used to extend the theory of dispersion to the case with thermal diffusion, and the result will depend on the choice of the closure for the heat-transfer problem.

The choice of Eqs ( 17) and ( 18), i.e., a local equilibrium theory, has already been explored in Davarzani et al. 2010(Davarzani et al., 2010a). We briefly give the main findings below. In the subsequent section, the original theory of thermal diffusion is presented in the case of a twoequation heat transfer model.

Mass Transport Model for Local Thermal Equilibrium

Application of the volume-averaging theorems on equation ( 5) and the use of classical simplifications lead to (Davarzani et al., 2010a)

( ) ( ) ( ) ( ) β β β β β β β β β βσ β β β β β β β β β β β βσ β β β β β β β β β β β ε ε ε ε ε ε ε βσ βσ v n n v . . . . . . . . 1 1 1 1 1 c T D dA T D T D c D dA c D c D c t c T A T T A ∇ - ∇ ∇ +         ∇ + ∇ ∇ + ∇ ∇ +         ∇ + ∇ ∇ = ∇ + ∂ ∂ - - - - - ∫ ∫ V V (30)
in which the spatial deviation of the concentration is defined by

β β β β c c c + = , in β V ( 31
)
The quasi-steady closure problem for the spatial-deviation concentration is obtained as

( ) ( ) β β β β β β β β β T D c D c c T . . . . ∇ ∇ + ∇ ∇ = ∇ + ∇ v v (32) BC1: ( ) ( ) β β β β β β βσ β β β β βσ T D c D T D c D T T ∇ + ∇ = ∇ + ∇ - . . n n , at βσ A (33) BC2: , i=1, 2, 3 (34)
Here, the point-concentration deviation is proportional to the macroscopic-average Problem IIa

β β β β β C C D b v b v 2 . ∇ = + ∇ (36)
BC:

β βσ β β βσ D D C n b n = ∇ - .
, at βσ A (37) Periodicity:

( ) ( ) r b r b β β C i C = +  , i=1,2,3 (38)
Averages:

0 = β β C b (39) Problem IIb β β β β β β T T S S D D b b b v 2 2 . ∇ + ∇ = ∇ (40)
BC:

( )

β βσ β β β β βσ T T T S D D D . . n b b n = ∇ + ∇ - , at βσ A (41) Periodicity: ( ) ( ) r b r b β β S i S = +  , i=1, 2, 3 (42)
Averages:

0 = β β S b (43)
We must note here that the closure problem for diffusion (Problem IIa) can be solved separately from the other closure problems. However, the closure problem for thermal diffusion The Darcy-scale mass balance equation with the thermal diffusion effect is then obtained as

(

) ( )

T c c t c T ∇ + ∇ ∇ = ∇ + ∂ ∂ . . . . D D v β β β β β β β β β β β β β β ε ε ε ε (44)
where the effective diffusion coefficients * β D and effective thermal diffusion coefficients * D β T are calculated from

β β β β βσ β β β βσ C A C dA V D b v b n I D 1 * -         + = ∫ (45) β β β β βσ β β β βσ β β β βσ βσ S A T T A S T dA V D dA V D b v b n I b n D * 1 1 -         + +         = ∫ ∫ (46)

Mass-Transport Model for Local Thermal Non-Equilibrium

Additional source terms in Eq. ( 35) must be taken into account in the case of local thermal nonequilibrium conditions. Following the approximations used before, we seek the concentration deviation under the form

( ) σ σ β β β σ σ βσ β β ββ β β β β T T r T T c c S S C - - ∇ + ∇ + ∇ = . . . ~b b b (47)
Here, we have supposed that the concentration deviation with respect to the averaged value is a function of concentration and temperature gradients in the fluid phase, the temperature gradient in the solid phase, and also the difference between the two averaged temperatures.

By substituting this new concentration deviation, and the β-phase temperature deviation from Eq. ( 25), into the quasi-steady closure problem presented before in Eqs ( 32)-(34), we obtain the closure problems for coupled heat and mass transfer listed as follows Problem IIIa

β β β β β C C D b v b v 2 . ∇ = + ∇ (48) BC: β βσ β β βσ D D C n b n = ∇ - .
, at βσ A (49) Periodicity:

( ) ( ) r b r b β β C i C = +  , i=1, 2, 3 (50)
Averages: 0 =

S i S = +  , i=1, 2, 3 (58)
Averages:

0 = β βσ S b (59)
Problem IIId

β β β β β β s D r D r T 2 2 . ∇ + ∇ = ∇ v (60)
BC:

( ) 0 . = ∇ + ∇ - β β β β βσ s D r D T n , at βσ A (61) Periodicity: ( ) ( ) r r β β r r i = +  , i=1, 2, 3 (62)
Averages:

0 = β β r ( 63
)
By substituting β c ~ and β T ~ from the decomposition given by Eqs ( 25) and (47) into Eq. ( 30), the Darcy-scale mass balance equation for the non-equilibrium two-equation temperature model case can be written

            ∇ + ∇ + ∇ ∇ =             - ∇ -       ∇ + ∂ ∂ σ σ βσ β β ββ β β β β σ σ β β β β β β β β β β β ε ε ε T T c T T c t c T T C . . . . . . * * * D D D u v (64)
where the effective tensors are defined by

β β β β βσ β β β βσ C A C dA V D b v b n I D 1 * -         + = ∫ (65) β ββ β ββ βσ β β ββ βσ β β ββ βσ βσ S A T T A S T dA V D dA V D b v b n I b n D * 1 1 -         + +         = ∫ ∫ (66) β βσ β βσ βσ β β βσ βσ β β βσ βσ βσ S A T T A S T dA V D dA V D b v b n b n D * 1 1 -         +         = ∫ ∫ (67) β β β β βσ β β β βσ β β β βσ βσ r dA s V D dA r V D A T A C v n n u 1 1 -         +         = ∫ ∫ (68)
We need therefore to solve four closure problems for concentration (Problem IIIa through IIId) after three closure problems for the two temperatures (Problem A-I through A-II) to model, at the Darcy-scale, mass transport coupled with the local non-equilibrium model defined in Eqs (64), ( 27) and (28).

The problem for β C b is the same as the problem IIa for the thermal equilibrium one-equation model, as it does not depend on the temperature gradient, which is consistent with the fact that we should recover the classical dispersion case if the Soret effect is zero.

In the next section, effective mass-transport coefficients are computed for simple unit cells in the case of the proposed two-equation model.

NUMERICAL RESULTS

Stokes equations are solved numerically over the unit cell shown in Fig. 2, representing the unit cell of a 2D periodic array of disks. The velocity field is generated by the introduction of an average pressure gradient along the x-axis. Given the velocity and deviation fields, the three closure problems for temperature can be solved independently from the mass-transport closure problems. As we have mentioned before, the closure problem for diffusion (problem IIIa) can also be solved independently. The closure problems have been solved numerically using COMSOL Multiphysics ® , a commercial package that is based on the finite element method. We used Lagrange quadratic elements for the velocity and linear elements for the pressure (P2-P1) to discretize the Stokes equation by considering Galerkin Least Square (GLS) streamline stabilization method [Hauke and Hughes, 1994]. The closure problems (in a form of diffusion or diffusion-advection equations) were solved using a quadratic Lagrange scheme. The simulation domains were discretized by triangular meshes with a minimum of 2,780 elements. Local mesh refinement was applied at the fluid-solid interface. To ensure that the output results are independent of the mesh size and structure, several mesh-sensitivity simulations were performed. In this study, all results are presented for a β-phase volume fraction, β ε , equal to 0.80. This value was chosen to be comparable to the results for the effective thermal diffusion coefficients with local-thermal equilibrium (Davarzani et al., 2010a).

The results are expressed as a function of the pore-scale thermal and solutal cell Péclet numbers, which are defined respectively as

𝑃𝑒 𝑇 = �𝜌𝑐 𝑝 � 𝛽 〈𝑣 𝛽 〉 𝛽 ℓ 𝑈𝑈 𝑘 𝛽 ( 69
)
𝑃𝑒 𝐶 = 〈𝑣 𝛽 〉 𝛽 ℓ 𝑈𝐶 𝐷 𝛽 (70)
where UC  is defined in Fig. 2.

Mass-Dispersion Tensors

The most important mass-transport tensor in Eq. ( 64) is the well-known dispersion tensor 𝐃 𝛽 * .

Problem IIIa has been solved numerically over the two-dimensional representative symmetric unit cell shown in Fig. 2 to give the mapping vector 𝐛 𝐶𝛽 . Then the dispersion tensors were calculated by using Equation ( 65), through the use of numerical integration on the fluid-solid interface and also the fluid subdomain. It is noteworthy to mention that, while the closure problems were solved in this section over a very simple symmetrical unit cell, the model introduced in this paper is valid for any kind of unit-cell geometry including non-symmetrical cells. Therefore, the proposed model can be applied, for example, to geometries coming from Xray tomography using a suitable method to correctly capture the anisotropy [START_REF] Guibert | A Comparison of Various Methods for the Numerical Evaluation of Porous Media Permeability Tensors from Pore-Scale Geometry[END_REF].

Our goal here was to understand the impact of the various pore-scale parameters on the resulting effective properties and not to provide a series of exhaustive results for various types of porous materials.

The longitudinal mass-dispersion coefficient, �𝐷 𝛽 * � 𝑥𝑥 , obtained for an averaged velocity along the x-axis, is shown in Fig. 3 Once again, a different behavior would be expected for a more disordered porous medium.

These dispersion curves can be represented as (Quintard and Whitaker, 1994) (while this formulation may lack accuracy in the transition regime)

�𝐷 𝛽 * � 𝑥𝑥 𝐷 𝛽 � = 1 𝜏 + 𝑎𝑃𝑒 𝐶 𝑛 (71)
where, the tortuosity term is defined as

1 𝜏 𝐈 = 𝐈 + 1 𝑉 𝛽 � 𝐧 𝛽𝜎 𝐛 𝐶𝛽 𝐴 𝛽𝛽 𝑑𝑑 ( 72
)
The exponent on the Péclet number, n, may vary between one for a linear-dispersive regime and a maximum value of two for a Taylor dispersion.

Thermal Diffusion Tensors

The closure problems IIIb and IIId have been solved over the same symmetric unit cell to obtain the thermal diffusion tensors 𝐃 𝑇𝛽𝛽 * and 𝐃 𝑇𝛽𝜎 * from the expressions presented in Eqs ( 66) and

(67).

The longitudinal component of the thermal diffusion tensor 𝐃 𝑇𝛽𝛽 * is plotted in Fig. 5 versus the Péclet number and for thermal conductivity ratio ranging from 0.01 to 100. To show the effective values in Fig. 5 for the diffusive regime more precisely, the results are subtracted from

1 𝜏 ⁄ = 1 1.20 ⁄ .
This value corresponds to the ratio of the effective thermal diffusion (or diffusion) coefficients to the (pore-scale) binary coefficients in the diffusive regime. Therefore, as shown in Fig. 5, this term is approximately zero for low Péclet numbers. When the Péclet number is increased, this term decreases and then it becomes negative for high Péclet numbers.

Similar behavior was also observed for the local-equilibrium case. (See thorough discussion in (Davarzani et al., 2010a)).

As can be seen in Fig. 5, the impact of the thermal conductivity ratio (𝜅 = 𝑘 𝜎 𝑘 𝛽 ⁄ ) on the effective thermal diffusion is important only in the convective regime. The effects are more important for a low-thermal conductivity ratio between 0.1 and 10. Outside that range, the effective thermal diffusion coefficients are less affected by variations. This is coherent with how thermal-dispersion tensors change with the thermal conductivity ratio as shown in Fig. A1a of Appendix A. In fact, the variation of the second term of the right-hand side of equation ( 65) with the thermal conductivity ratio influences the effective thermal diffusion coefficient. We should note again that the two-equation model for temperature may not predict the temperature correctly when the solid/fluid conductivity ratio is very low and the solid diameter of porous media is large. In this case, we observe a difference between solid center and solid surface temperature [START_REF] Debenest | Transport in Highly Heterogeneous Porous Media: From Direct Simulation to Macro-Scale Two-Equation Models or Mixed Models[END_REF][START_REF] Xin | Numerical Investigation of Local Thermal Non-Equilibrium Effects in Coal Porous Media with Cryogenic Nitrogen Injection[END_REF]. Therefore, a mixed model can also be used if more accuracy is required and if there is no restriction in the computational cost.

The transverse component of the thermal diffusion tensor 𝐃 𝑇𝛽𝛽 * is represented in Fig. 6 versus the Péclet number and also for different values of the thermal conductivity ratio. The second thermal diffusion tensor appearing in the macroscopic-mass transport, Eq. ( 64), is , the impact of the thermal conductivity ratio is greater for moderate and high thermal conductivity ratios (the change between the thermal conductivities of 10 and 100 is considerable in Fig. 7).

𝐷 𝑇𝛽𝜎 * .
The transverse component of the thermal diffusion tensor �𝐷 𝑇𝛽𝜎 * � 𝑦𝑦 is plotted in Fig. 8 as a function of the Péclet number and for different thermal conductivity ratios. We see that this

for all ranges of the Péclet number and thermal conductivity ratio.

Mass-Convective Vector

Simulations of Problem IIId have been obtained to calculate the vector 𝐮 𝐶𝛽 in Eq. ( 64), which was named "the mass-convective vector" to distinguish it from "the heat-convective vectors" presented in Appendix A for the heat-transfer non-equilibrium model. The solution of Problem

IIId for 𝑟 𝛽 depends on the mapping vector 𝑠 𝛽 of the thermal non-equilibrium closure problems.

Therefore, Problem A-III must be solved before or simultaneously with Problem IIId to obtain the mapping vector 𝑟 𝛽 . Then the mass-convective vector is calculated using the formula given in Eq. ( 68). for moderate thermal conductivity ratios. However, increasing the thermal conductivity ratios has less effect on mass-transport-convective term for high thermal conductivity ratios.

Soret Coefficients

To evaluate the importance of the effective thermal diffusion tensors, compared to dispersion, we rewrite Eq. ( 64) in term of Soret coefficients as

( ) ( ) ( ) ( ) ( ) σ σ βσ β β ββ β β β β σ σ β β β β β β β β β β β ε ε ε T T c T T c t c T T C ∇ + ∇ + ∇ ∇ = - ∇ - ∇ + ∂ ∂ . . . .

S S

). However, this is no longer valid when convection becomes important. In that case, the Soret ratio decreases gradually with increasing Péclet number and then becomes negative with asymptotic behavior, for the convection-dominant regime. The impact of the thermal conductivity ratio is more important for moderate thermal conductivity ratios between 0.1 and 10. Similar to longitudinal-mass dispersion, the transition between convective and diffusive regime occurs when the Péclet number is around one.

The results for the transverse-Soret coefficient, ( ) Finally, the relationship between the effective properties obtained from non-equilibrium and equilibrium models can be also expressed in term of the thermal-equilibrium Soret coefficients, 𝐒 𝑇𝛽 * , and non-equilibrium Soret coefficients as

𝐒 𝑇𝛽 * = 𝐒 𝑇𝛽𝛽 * + 𝐒 𝑇𝛽𝜎 * (75)

CONCLUSIONS

We have developed a Darcy-scale theoretical model to predict coupled heat and mass transfer in porous media including the Soret cross-effect, for local thermal non-equilibrium regimes.

The porous-medium energy and mass-transfer equations were deduced using the volumeaveraging upscaling method. The model equation for mass balance contains a dispersion tensor, two thermal diffusion tensors involving the gradient in the fluid and solid phase, and also a mass-convective vector corresponding to the temperature gradient of the temperature difference between the two phases.

Effective Darcy-scale coefficients have been determined by solving the corresponding closure problems over a simple representative-unit cell. The results show that, except for pure diffusive cases where the Soret coefficients in the porous medium and the free fluid are the same, the effective-Soret coefficients are impacted by the fluid hydrodynamics and also depend significantly on the medium's thermo-physical properties, in particular, the phases' thermal conductivity contrasts.

When the conditions of local-thermal equilibrium are not satisfied, it may be more convenient and accurate to use the theoretical model proposed in this study for modeling coupled transfer in porous media with a thermal diffusion-cross effect, especially at large Péclet numbers. The proposed model is original and must be tested thoroughly; in particular, experimental studies must be developed to validate the theoretical model.

closure problems for the thermal dispersion case. The three closure problems are listed below.

Problem A-I ( ) ( ) ββ β ββ β β β ββ β β ε ρ ρ c b v b v 1 2 . - - ∇ = + ∇ T p T p k c c , in β V (A. 1) BC1: β βσ σβ σ βσ ββ β βσ k k k T T n b n b n - ∇ = ∇ . . , at βσ A (A. 2) BC2: σβ ββ T T b b = , at βσ A (A. 3) σβ σ σβ σ ε c b 1 2 0 - - ∇ = T k , in σ V (A. 4) Periodicity: ( ) ( ) r b r b ββ ββ T i T = +  & ( ) ( ) r b r b σβ σβ T i T = +  , i=1, 2, 3 (A. 5) Averages: 0 = ββ T b & 0 = σβ T b (A. 6) dA k A T ∫ ∇ = βσ ββ β βσ ββ b n c . 1 V (A. 7) ββ σβ σ σβ σβ βσ c b n c - = ∇ = ∫ dA k A T . 1 V (A. 8) Problem A-II ( ) βσ β βσ β βσ β β ε ρ c b b v 1 2 . - - ∇ = ∇ T T p k c
, in β V (A. 9) BC1:

σ βσ σσ σ βσ βσ β βσ k k k T T n b n b n + ∇ = ∇ . . , at βσ A (A. 10) BC2: σσ βσ T T b b = , at βσ A (A. 11) σσ σ σσ σ ε c b 1 2 0 - - ∇ = T k , in σ V (A. 12) Periodicity: ( ) ( ) r b r b βσ βσ T i T = +  & ( ) ( ) r b r b σσ σσ T i T = +  , i=1, 2, 3 (A. 13) Averages: 0 = βσ T b & 0 = σσ T b (A. 14) dA k A T ∫ ∇ = βσ βσ β βσ βσ b n c . 1 V (A. 15) βσ σσ σ σβ σσ βσ c b n c - = ∇ = ∫ dA k A T . 1 V (A. 16) Problem A-III ( ) β β β β β β β ε ρ h s k s c p 1 2 . - - ∇ = ∇ v , in β V (A. 17) BC1: σ σ βσ β β βσ s k s k ∇ = ∇ . . n n
, at βσ A (A. 18) BC2:

σ β s s + = 1 , at βσ A (A. 19) σ σ σ σ ε h s k 1 2 0 - + ∇ = , in σ V (A. 20) Periodicity: ( ) ( ) r r β β s s i = +  & ( ) ( ) r r σ σ s s i = +  , i=1, 2, 3 (A. 21) Averages: 0 = β s & 0 = σ s (A. 22) dA s k h A ∫ ∇ = βσ β β βσ β . 1 n V (A. 23) T A h h dA s k h = = ∇ - = ∫ β σ σ σβ σ βσ . 1 n V (A. 24)
We note that the volumetric heat-transfer coefficient, T h , can be calculated from solving Problem A-III and using equation (A. 23) or (A. 24). The other effective coefficients are given below.

( ) The results of this model for ɛ 𝛽 = 0.38 (results have not been presented in this study) and thermal-dispersion coefficients agree quite well with the results obtained in [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF]. Generally, the influence of the phases' thermal conductivity contrasts on the effective thermal-dispersion coefficients is greater for moderate-thermal conductivity ratios between 0.1 and 10. As mentioned before, this behavior can be related to corresponding thermal diffusion coefficients.

ββ β β ββ βσ β β ββ ρ ε βσ T p A T c dA k b v b n I k 1 -         + = ∫ V (A.
To calculate the convective heat vectors from Eqs (A. 29) to (A. 32) we must solve all three closure problems for thermal dispersion (A-I, A-II, and A-III). The results are shown in Fig. A2 as a function of Péclet number and for ɛ 𝛽 = 0.80 and different thermal conductivity ratios.

Except for the high Péclet numbers, the results agree with the results presented in [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF] for ɛ 𝜷 = 0.38. Their results show an absolute maximum value around Pe=300 followed by asymptotic behavior toward zero for Péclet numbers approaching infinity. By contrast, the convective vectors continue to increase in this study for both porosities (ɛ 𝜷 = 0.38 and 0.80)

without reaching an absolute maximum value. Below this Péclet number, the results agree with the results obtained in Quintard et al. (1997) † .

The results for the heat transfer coefficients are plotted in Fig. A3 as a function of Péclet number and for different thermal conductivity ratios are in good agreement with the ones published before in [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF]. thermal conductivity ratios, 𝜅, and 𝜀 𝛽 = 0.80 † A minor error in the Fortran code (M. Quintard, personal communication) used to calculate the integrals giving these coefficients may explain these differences and suggests that these results should be considered instead of those in [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF].

KEY WORDS :

 : Coupled phenomena, Soret effect, Numerical simulation, Upscaling, Volume averaging NOMENCLATURE Area of the β-σ interface contained within the macroscopic regionaverage mass fraction in the β-phase T Pe Thermal cell Péclet number Spatial deviation mass fraction in the β-phase C Pe Solutal cell Péclet number Binary diffusion coefficient, m 2 /s Radius of the averaging volume, m Thermal diffusion coefficient, m 2 /s.K Mass average velocity in the β-phase, m/s Thermal conductivity of the solid phase, β ε Volume fraction of the β-phase or porosity phase equation, W/m.K Total mass density in the β-phase, kg/m 3

T

  are the temperature fields in the β-and σ-phases, respectively, β k and σ k are the fluid and solid thermal conductivity coefficients, βσ n is a unit normal vector directed from the fluid toward the solid phase, and βσ A is the area of the interface between the two phases (β-σ interfaces).

b

  are the closure variables for temperature.

  Because of scale limitations (time and length), the macroscopic behaviors of porous media cannot generally be represented by a one-equation model for two different phases. The twoequation model consists of separate heat-transport equations for both the σ-and β-phases. The dominant coupling between the two-phase energy equations is then represented by a flux between the phases that depends, in the simpler version, on an exchange coefficient and the difference between the temperatures of the two phases. The temperature spatial deviations in terms of the macroscopic source terms are approximated in this case as scalars. They can be obtained from the solution of pore-scale closure problems as detailed in

  variables are the solution of the two boundary value problems IIa and IIb presented below.

  ) depends on the closure variable for temperature β T b .

  The three remaining closure problems depend on the closure problems for temperature. Here, the solution of the closure problem IIIb depends on the solution of the closure problem A-I for 𝐛 𝑇𝛽𝛽 ; the closure problem IIIc depends on 𝐛 𝑇𝛽𝜎 , which comes from the solution of the closure problem A-II; and, finally, the solution of the closure problem IIId depends on the solution of the closure problem A-III for the scalar value, 𝑠 𝛽 .

  as a function of the Péclet number. The dispersion tensor features the usual behavior with three different regimes: a diffusion-dominant regime for low Péclet numbers, a transition behavior for moderate Péclet numbers, and finally, a convection-dominant regime with a power-law trend. Changes in the transverse component of the mass dispersion coefficient, �𝐷 𝛽 * � 𝑦𝑦 , with the Péclet number is shown in Fig. 4. Given the symmetries of the unit cell, the transverse coefficient is equal to the longitudinal coefficient in the diffusive-dominant regime i.e., the medium is indeed isotropic. The transition regime for the transverse component occurs at a Péclet number greater than for the longitudinal component. For such a simple unit cell, the longitudinal components are much bigger than transverse components in the convective regime.

  For this symmetric unit cell, the transverse component of the thermal diffusion coefficient �𝐷 𝑇𝛽𝛽 * � 𝑦𝑦 is constant and equal to the longitudinal thermal diffusion coefficient, �𝐷 𝑇𝛽𝛽 * � 𝑥𝑥 , in the diffusive regime. As expected, the results show that the transverse thermal diffusion coefficient is smaller than the longitudinal coefficient in the convective regime. Here also, like for the thermaldispersion tensors, the variations of the thermal conductivity ratio change more considerably than the effective coefficients for moderate thermal conductivity ratios. The transverse thermal diffusion coefficient, �𝐷 𝑇𝛽𝛽 * � 𝑦𝑦 , decreases with increasing Péclet number in the convective regime.

  The longitudinal component of this thermal diffusion tensor, �𝐷 𝑇𝛽𝜎 * � 𝑥𝑥 , is shown in Fig.7versus the Péclet number and for different thermal conductivity ratios. It is clear that for the diffusion-dominant regime, this effective thermal diffusion value is close to zero. Increasing the Péclet number increases the effective coefficient: lightly for the low thermal conductivity ratios and significantly for high thermal conductivity ratios. These results are also coherent, in terms of behavior, with the results for the effective thermal conductivity coefficient,

  component also is equal to its longitudinal coefficient for low Péclet number, i.e., zero values for the diffusion-dominant regime, and then that the effective values increase in the convectiondominant regime for the moderate and high thermal conductivity ratios. The effective values for the low contrast of phase-thermal conductivities are much less sensitive to changing Péclet number. The transverse effective values are lower than the longitudinal values, and the transition between diffusion-and convection-dominant regime occurs when the Péclet number is around 10. Here also, the behavior results are comparable with �𝑘 𝛽𝜎 * � 𝑦𝑦 as shown in Fig. A1d of Appendix A. We found from the numerical results that the thermal diffusion tensors come from the local thermal non-equilibrium, 𝐃 𝑇𝛽𝛽 * and 𝐃 𝑇𝛽𝜎 * , and the thermal diffusion tensors obtained from the local thermal equilibrium, 𝐃 𝑇𝛽 * , are related by 𝐃 𝑇𝛽 * = 𝐃 𝑇𝛽𝛽 * + 𝐃 𝑇𝛽𝜎 *

  The x-component of vector 𝐮 𝐶𝛽 is plotted in Fig.9. In this figure, the value �𝑢 𝐶𝛽 � 𝑥 has been divided by the average velocity 〈𝑣 𝛽 〉 𝛽 , to evaluate the significance of the mass-convective term compared to the mass-transport-convective term. The mass-convective term versus Péclet number gives a reverse sigmoid curve. The results show a constant and negligible masstransport-convective term for a low thermal conductivity ratio. The mass-conductive terms are constant in the diffusive regime and then decrease in the conductive regime. They are significant

  Fig. 10 and Fig. 11 represent the variations of these two Soret coefficients versus the Péclet

  Fig.10b, are quite similar.Here also, the effective-transverse-Soret coefficient in porous media ( ) the free fluid, in the purely diffusive regime ( ( ) coefficient decreases with increasing Péclet number. The transition regime occurs when the Péclet number is around 1. No asymptotic behavior has been observed for high Péclet numbers, instead absolute minimum values were found for low thermal conductivity ratios. The effective-transverse-Soret coefficient ( )

Fig. 11

 11 Fig. 11 shows the ratio of the effective-Soret coefficients, 𝐒 𝑇𝛽𝜎 * , versus the Soret coefficient in

  developed using COMSOL Multiphysics ® , following the same numerical procedures outlined in Section 4. Longitudinal and transverse thermal-dispersion coefficients are plotted as a function of thermal Péclet number for ɛ 𝛽 = 0.80 and for different thermal conductivity ratios in Fig. A1. They are calculated from the numerical solution of Problem A-I and A-II using Eqs (A. 25)-(A. 28).

FIGFIG. 5 :FIG. 6 :FIG. 7 :

 567 FIG. 1: Problem configuration

FIG. 8 :

 8 FIG. 8: Transverse thermal diffusion coefficient �𝐃 𝐓𝛃𝐓 * � 𝐲𝐲 versus Péclet number for different

APPENDIX A. THE CLOSURE PROBLEMS FOR THE LOCAL THERMAL NON-EQUILIBRIUM MODEL: EFFECTIVE CALCULATION RESULTS AND MODEL

VALIDATION

The closure problems for the two-temperature model have been already investigated [START_REF] Carbonell | Heat and Mass Transfer in Porous Media, Fundamentals of Transport Phenomena in Porous Media[END_REF][START_REF] Levec | Longitudinal and Lateral Thermal Dispersion in Packed Beds. Part I: Theory[END_REF][START_REF] Quintard | Diffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF][START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF][START_REF] Vortmeyer | Equivalence of One-and Two-Phase Models for Heat Transfer Processes in Packed Beds: One Dimensional Theory[END_REF][START_REF] Zanotti | Development of Transport Equations for Multiphase System-I: General Development for Two Phase System[END_REF]. Here, we present these three closure problems and the corresponding macroscopic-effective coefficients needed to solve the Darcy-scale heat-transfer equations presented in Eqs ( 27) and (28).

The numerical results are obtained for the same unit cell as in Fig. 2. We followed the decomposition method proposed by [START_REF] Quintard | Two-Medium Treatment of Heat Transfer in Porous Media: Numerical Results for Effective Properties[END_REF] to solve the integro-differential


Daniel Chatroux

Julien Chauvin

Utilisation des MOSFET en mode avalanche

Keywords: Avalanche, écrêtage, MOSFET, mise en série

RESUME

-Le phénomène d'avalanche des MOSFET permet d'écrêter la tension par le composant lui-même, plutôt que de lui associer une diode dite Zéner ou une fonction écrêteur par le circuit. La présentation est centrée sur l'évolution des spécifications, des méthodes de caractérisation et des performances des composants en avalanche. Une amélioration d'une méthode de caractérisation est proposée. D'autre part l'utilisation de l'avalanche dans des applications est illustrée par un commutateur à MOSFET 100 kV utilisé depuis vingt ans sur le Tokamak Tore Supra à Cadarache et par le développement récent d'un commutateur 10 kA pour le test en court-circuit à durée contrôlée des accumulateurs ou batteries Lithium ion.

INTRODUCTION

Très souvent, l'utilisation de l'avalanche des transistors pour écrêter une tension est considérée comme dangereuse, préjudiciable à la durée de vie des composants, et non utilisable pour des transistors en parallèles. Dans les années 2000, le CEA a démontré qu'il était possible d'utiliser des MOSFET caractérisés en courant d'avalanche et montés en parallèle et en série sur des commutateurs de hautes tensions et forts courants avec d'excellentes performances. L'objectif de cette présentation est de faire un état des lieux des connaissances d'aujourd'hui sur l'avalanche, de l'évolution des spécifications et des méthodes de mesures aisément disponibles pour les concepteurs.

AVALANCHE DES MOSFET

Phénomène physique de l'avalanche des composants

L'avalanche correspond au passage de courant lors de l'atteinte d'une tension correspondant au champ électrique critique dans la jonction tenant la tension, avec maintien de la tension aux bornes.

C'est ce phénomène qui est utilisé dans les diodes dites Zéner, ainsi que les diodes écréteuses de puissance comme les Transils®. L'effet Zener est un écrêtage à un seuil de quelques volts seulement. Au-delà, l'écrêtage s'effectue par l'effet d'avalanche. Pour des tensions voisines de 5-6V ces effets à coefficients de température opposés coexistent et permettent une diminution de la dérive en température qui est connue. Pour les diodes écrêteuse de puissance, il n'y a pas de spécification de niveau maximum de courant d'avalanche. Par contre, les contraintes thermiques imposent une énergie maximale pour les durées courtes et une puissance maximale que peut dissiper le composant pour les durées longues.

La physique de l'avalanche des MOSFET est aujourd'hui bien décrite dans certaines notes d'applications [START_REF]Power MOSFET Avalanche Design Guidelines[END_REF]. Pour les MOSFET, la diode interne est en fait un transistor bipolaire interne dont la base et l'émetteur sont court-circuités. En mode avalanche, il existe un risque qui est d'amorcer localement le transistor par le phénomène de second claquage, ce qui conduirait à sa destruction par point de fusion localisé. C'est donc le dépassement d'un niveau de courant qui donne lieu à la destruction, selon la qualité de maintien du court-circuit sur toute la surface du composant.

Evolutions des spécifications au cours du temps

Les MOSFET n'étaient classiquement pas spécifiés en mode avalanche, ou spécifiés en énergie d'avalanche pour des niveaux de courants faibles, bien inférieurs aux besoins des applications. Il y a vingt ans, les travaux du CEA sur les commutateurs de haute tension ont montré que le critère dimensionnant était le niveau de courant [START_REF] Chatroux | Power switch : the standard small components strategy[END_REF]. Par exemple, pour des durées d'avalanche de 200 ns et 1µs, les composants cassaient exactement au même niveau de courant d'avalanche. Cela signifie bien que l'énergie (qui variait alors d'un facteur 5) n'est pas le critère pertinent. C'est le niveau de courant qui est le paramètre limitant. D'autre part, c'est le design du composant qui assure un bon court-circuit base émetteur du transistor parasite du MOSFET. Les travaux visaient à sélectionner et qualifier une référence de composant d'un fournisseur, puis à utiliser les composants sans aucun tri individuel. Il fallait par contre surveiller que le design n'évolue pas chez le fabriquant au cours du temps.

Ceci avait permis de tracer l'aire de sécurité avec l'avalanche de l'IRF840 de IR, qui est un MOSFET 8A 500V. Le composant supportait 60A en avalanche. L'aire de sécurité est très étendue, mais il restait encore une zone où il peut être détruit, par exemple pour l'ouverture d'un courant de courtcircuit supérieur à 60A. Une limitation de la tension de grille à 10 volts permettait de rester dans cette aire.

Aujourd'hui, les spécifications des MOSFET en avalanche ont beaucoup évoluées. Certains composants sont spécifiés en mono-coup, mais aussi en mode répétitif pour un courant d'avalanche semblable à leur courant nominal. Lors de la conception, il faut alors seulement ajouter les pertes en phase d'avalanche, lors du calcul des pertes de commutation.

Caractérisation des composants en mode avalanche

Historiquement au CEA, la caractérisation en avalanche des MOSFET se faisait à courant constant. Ce courant carré facilitait les interprétations. Par contre le banc de test (Figure 2) nécessite une alimentation et un interrupteur électronique de haute tension associés à une résistance non inductive de forte énergie. Si la mesure de courant ne pose pas de problème de métrologie (sondes à effet Hall sur oscilloscopes), une mesure de tension fidèle, précise et rapide aux bornes du MOSFET est délicate à partir de quelques centaines de volts. La solution proposée est d'utiliser un diviseur entre une résistance 4950 ohms (constituée de résistances couche métalliques en séries pour tenir la tension à l'ouverture), placées en série avec un câble 50 ohms chargé par 50 ohms sur l'entrée de l'oscilloscope (Figure 4). Les mesures basses impédances, de l'ordre de 5 kiloohms, sont connues dans le domaine des hautes tensions pulsées lorsqu'on veut étudier des impulsions ou des fronts rapides. Cette technique est ici utilisable du fait de la faible tension d'alimentation. Le concepteur peut alors étudier aisément le comportement en avalanche des composants, rechercher les limites en courants et pour certains composants observer leur tendance à osciller spontanément, qui induit un risque fort de casse en utilisation. Des références de diodes ultrarapides ont été caractérisées en avalanche. Le critère choisi a priori pour l'utilisation en série et matrice était que la diode supporte son courant nominal en avalanche. L'expérience a montré par la suite la pertinence de ce choix.

Le commutateur 25kV 2kA utilisait deux blocs en parallèles constitués de cinq cartes 5kV 1kA, mises en séries. Cette technologie a été réutilisée pour les interrupteurs de sécurité du Tokamak Tore Supra à Cadarache, rebaptisé WEST. En cas de court-circuit dans l'un des Klystrons qui fournissent des microondes pour chauffer le plasma, l'interrupteur 100 kV 20A ouvre le courant de défaut [START_REF] Beaumont | 100 kV solid-state switch for fusion heating systems[END_REF]. Vingt ans plus tard, certains des équipements développés sont encore utilisés, comme celui présenté ci-dessus (Figure 8).

Commutateur basse tension 10 kA pour test de

batteries Lithium-ion Le CEA mène des recherches sur les batteries Lithium-ion depuis les matériaux jusqu'aux applications transport, c'est-àdire la conception des batteries répondant aux besoins des véhicules électriques routiers, maritimes ou aériens. Pour l'étude de la sécurité de ces batteries, qui est un des points durs de la technologie, le CEA s'est muni de plates-formes de tests abusifs. Pour ces besoins de tests, un commutateur de mise en courtcircuit a été développé pour un courant maximum de 10 kA.

Pour cette application, du fait du niveau de courant, le principal problème à traiter est l'énergie des quelques mètres de câblage entre la batterie en test et l'interrupteur électronique. Le risque de feu impose une distance minimale de 1,50 mètre. Malgré un câblage conçu pour être faiblement inductif, l'énergie des câblages à absorber est de 50 Joules. La résistance équivalente doit être inférieure au milli-ohm.

Plutôt que d'utiliser des composants de plus haute tension associé à un écrêteur, il a été choisi d'optimiser le système en utilisant 100 interrupteurs MOSFET 100 V en parallèle et d'écrêter l'énergie dans l'ensemble de ces composants.

Une première phase a permis de vérifier la tenue en courant d'avalanche de trois MOSFET en parallèle. Ensuite une carte 2 kA a été réalisée : A cette occasion, il a été montré que le paramètre pertinent était le courant avec un seuil à ne pas dépasser. La limite en énergie est la conséquence thermique selon la durée de test choisie.

Pour les diodes, il avait été choisi arbitrairement de sélectionner les références de diodes supportant leur courant nominal en avalanche pour l'utilisation en série et en matrice. Ce choix s'est avéré pertinent dans les années qui ont suivi.

Aujourd'hui, les MOSFET ne sont toujours pas spécifiés en courant maximum d'avalanche, mais sont spécifiés en énergie pour des niveaux de courant élevés. Ce niveau renseigne sur la robustesse du composant en avalanche.

Avec les composants spécifiés 100A en avalanche, un interrupteur 10 kA apte à absorber l'énergie des câblages avec un module Lithium à tester a été réalisé.

Figure 1 :

 1 Figure 1: Aire de sécurité de l'IRF840

Figure 2 :

 2 Figure 2 : Banc de test en avalanche à courant constant Pour les spécifications actuelles des MOSFET, le montage de test utilise l'ouverture d'un courant d'un circuit inductif [1], Figure 3.

Figure 3 :

 3 Figure 3: Test d'un MOSFET en avalanche avec forme de courant triangulaire.

Figure 4 :

 4 Figure 4 : sonde 5 kohm sortie 50 ohmsDu fait de la sortie en 50 ohms, la bande passante de la sonde réalisée peut être caractérisée en utilisant la sortie de calibration d'un analyseur de spectre (Figure5).

Figure 5 :

 5 Figure 5 : comportement fréquenciel de la sonde réalisée 3. EXEMPLES DE REALISATION 3.1. Interrupteur haute tension 100 kV à base de MOSFET en série Pour les besoins d'augmenter la fiabilité et de baisser le coût de fonctionnement de lasers impulsionnels devant fonctionner en continu pendant des années, le CEA a développé, il y a vingt ans des commutateurs électroniques de haute tension (25 kV) fort courant (2kA) commutant en 20 ns utilisant 3500 MOSFET en série et parallèle [2]. Du fait d'une commande de fort courant permettant de commander les composants d'un étage malgré un MOSFET en défaut, la technologie développée est tolérante à la panne et permet une durée de vie de 30 ans, sans besoin d'intervention [3]. Cette technologie utilise l'avalanche de MOSFETs qui sont non spécifiés pour cet usage, mais dont la référence a été caractérisée en courant d'avalanche, sans tri individuel.

Figure 6 :

 6 Figure 6 : Carte 5 kV 1kA Les cartes étaient placées au sein d'une cuve à huile minérale avec une circulation d'huile et échangeur externe huile/eau. Le commutateur était installé en lieu et place des commutateurs à gaz (thyratron) utilisés précédemment. Ces tubes avaient une durée de vie limitée à 1000 heures.

Figure 7 :

 7 Figure 7 : commutateur 25 kV 2 kA ouvert

Figure 8 :

 8 Figure 8 : Ensemble des Klystrons et des leurs interrupteurs associés

Figure 9 :

 9 Figure 9 : Carte 2 kA constituée de 20 MOSFET Les fusibles en série avec chaque MOSFET permettent d'ouvrir son courant en cas de défaillance de celui-ci.

Figure 10 :

 10 Figure 10 : test à 10 kA La figure 10 présente le test d'ouverture à 10 kA avec avalanche à la tension de 122V dans l'ensemble des MOSFET. Pour cette réalisation, les composants sont utilisés en respectant leur spécification de courant d'avalanche de 100 A. Pour des

Figure 11 :

 11 Figure 11 : courant d'avalanche maximum Ce MOSFET a un courant maximum en avalanche de 630A, qui est 5 fois son courant nominal ; 1,3 fois le courant pulsé maximum spécifié, et 6 fois le courant de test en avalanche.

Figure 12 :

 12 Figure 12 : tension d'écrêtage en fonction de la températue La résistance série du phénomène d'avalanche varie peu, elle est de 50 milliohms pour une température de 25°C et atteint 80 milliohms pour une température de 150°C. Cette valeur est à comparer à la RDSON du composant qui est de 2,3 milliohms. La forte résistance série apparente du phénomène d'avalanche tend à équilibrer les courants d'avalanche entre les MOSFET en parallèle. Pour des composants de plus haute tension, on a pu observer un courant d'avalanche plus limité pour un composant de type CoolMOS, et une tendance à osciller.

Figure 13 :Figure 14 :

 1314 Figure 13 : courant d'avalanche d'un CoolMOS 20A On a aussi pu observer le très bon comportement d'un MOSFET SiC.

  Au niveau des tests, comme le courant maximum d'avalanche est le paramètre pertinent, un test d'ouverture inductif est suffisant et très simple de mise en oeuvre. Une sonde basse impédance 5 kohms a été proposée pour la mesure précise de la tension L'évolution en température et la résistance apparente équivalente du composant choisi ont été présentés. Parmi les autres MOSFET de plus haute tension testés, les performances en avalanche d'un premier MOSFET SiC ont été caractérisées démontrant ainsi un bon comportement en avalanche de cette technologie.

Dans le cas de MOSFET en parallèle, comme dans notre réalisation, l'équilibrage des courant est défini par la dispersion de la tension d'écrêtage des composants, la résistance série apparente du phénomène d'avalanche et l'évolution de la tension d'avalanche avec la température. Nous avons mené des essais pour mesurer ces paramètres.