Francisco J Valdés-Parada

Didier Lasseux
email: didier.lasseux@ensam.eu

Fabien Bellet

A new formulation of the dispersion tensor in homogeneous porous media

Keywords: Dispersion, e↵ective di↵usion, volume averaging, closure problem

Dispersion is the result of two mass transport processes, namely molecular di↵usion, which is a pure mixing e↵ect and hydrodynamic dispersion, which combines mixing and spreading. The identification of each contribution is crucial and is often misinterpreted. Traditionally, under a volume averaging framework, a single closure problem is solved and the resulting fields are substituted into di↵usive and dispersive filters. However the di↵usive filter (that leads to the e↵ective di↵usivity) allows passing information from convection, which leads to an incorrect definition of the e↵ective medium coe cients composing the total dispersion tensor.

In this work, we revisit the definitions of the e↵ective di↵usivity and hydrodynamic dispersion tensors using the method of volume averaging. Our analysis shows that, in the context of laminar flow with or without inertial e↵ects, two closure problems need to be computed in order to correctly define the corresponding e↵ective medium coe cients. The first closure problem is associated to momentum transport and needs to be solved for a prescribed Reynolds number and flow orientation. The second closure problem is related to mass transport and it is solved first with a zero Péclet number and second with the required Péclet number and flow orientation. All the closure problems are written using closure variables only as required by the upscaling method. The total dispersion tensor is shown to depend on the microstucture, macroscopic flow angles, the cell (or pore) Péclet number and the cell (or pore) Reynolds number. It is non-symmetric in the general case. The condition for quasi-symmetry is highlighted. The functionality of the longitudinal and transverse components of this tensor with the flow angle is investigated for a 2D model porous structure obtaining consistent results with previous studies.

Pore-scale model

Dispersion of a solute (species A) in porous media is a fundamental subject that has been largely studied over the past century and it remains as an interesting study field due to the wide range of applications that it encompasses. This transport mechanism is the result of di↵usion, which is related to a pure mixing e↵ect, and variations in the convective fluxes within the pores leading to both mixing and spreading. The governing pore-scale mass balance equation for a solute (species A) transported within the fluid phase (the -phase) that saturates the porous medium is

@c A @t + r • (c A v ) = r • (D rc A ) , in the -phase (1)
Here c A and D respectively denote the pointwise species A concentration and the molecular di↵usivity, the latter being considered as a constant in this work, whereas v is the fluid velocity, which satisfies the total mass and momentum conservation equations at the pore-scale:

r• v = 0, in the -phase (2a) ⇢ v • rv = rp + µ r 2 v , in the -phase (2b)
Note that flow has been assumed to be steady, incompressible and Newtonian. Without loosing generality, gravity has been omitted in the momentum equation. In Eq. (2b), ⇢ and µ denote the fluid density and dynamic viscosity, respectively, which are assumed to be constants. Furthermore, for the sake of simplicity in the analysis, the porous medium is assumed to be rigid and homogeneous, so that intrinsic average properties (e.g. the porosity, the permeability, etc.) are position-invariant. In addition, the non-slip boundary condition is imposed at the solid-fluid interface

A v = 0, at A (3)
and the solid phase (i.e., the -phase) is assumed impermeable to mass transfer, so that

n • D rc A = 0, at A (4)
Despite the simple form of the governing equations at the pore-scale, all the essential elements are contained to give rise to the well-known convection-dispersion equation for mass transport after an upscaling process is applied. The study of dispersion from di↵erent theoretical points of view is available in several references, a couple of examples being the works by Cushman et al. [START_REF] Cushman | A primer on upscaling tools for porous media[END_REF] and Chapter 11 of Sahimi [START_REF] Sahimi | Flow and Transport in Porous Media and Fractured Rock, 2nd Edition[END_REF]. In particular, the study of dispersion using the volume averaging method dates back to the classic paper of Whitaker in 1967 [START_REF] Whitaker | Di↵usion and dispersion in porous media[END_REF], in which the macroscopic model was derived from the pore-scale equations by means of an averaging process. However, no closure procedure was provided at that moment and thus the e↵ectivemedium coe cients involved in the model had to be estimated experimentally. The model derivation was later improved by correctly defining the spatial decomposition of pore-scale variables as detailed by Gray [START_REF] Gray | A derivation of the equations for multiphase transport[END_REF]. The closure process, under a volume averaging approach, was presented in a series of works dealing with dispersion in pulsed systems. Firstly, transport of a solute in capillaries undergoing heterogeneous reaction and adsorption was presented by Paine et al. [START_REF] Paine | Dispersion in pulsed systems I. Heterogeneous reaction and reversible adsorption in capillary tubes[END_REF]. Their study showed the time dependence of the dispersion coe cient under di↵erent reactive conditions. As expected, under quasi-steady conditions, the well-known

Taylor [START_REF] Taylor | Dispersion of soluble matter in solvent flowing slowly through a tube[END_REF]-Aris [START_REF] Aris | On the dispersion of a solute in a fluid flowing through a tube[END_REF] result was recovered. Following this approach, the closure problem for dispersion in porous media was derived by Carbonell and Whitaker [START_REF] Carbonell | Dispersion in pulsed systems II: Theoretical developments for passive dispersion in porous media[END_REF] and it was subsequently solved in periodic unit cells by Eidsath et al. [START_REF] Eidsath | Dispersion in pulsed systems III: Comparison between theory and experiments for packed beds[END_REF] obtaining good agreement with experimental data for particle Péclet number values below 1000 [see Fig. 13 in 9]. This type of analysis was soon extended to study heterogeneous porous materials by Plumb and Whitaker [START_REF] Plumb | Dispersion in heterogeneous porous media: 1. Local volume averaging and largescale averaging[END_REF] and the corresponding closure problem was solved in spatially periodic unit cells for stratified porous media [START_REF] Plumb | Dispersion in heterogeneous porous media: 2. Predictions for stratified and two-dimensional spatially periodic systems[END_REF]. In a later work by Quintard et al. [START_REF] Quintard | Dispersion in heterogeneous porous media: One-equation nonequilibrium model[END_REF], a one-equation model, without involving large-scale mass equilibrium, was derived for studying mass dispersion in heterogeneous porous media. In addition, solutions to the closure problems have also been carried out for two-equation models in fissured media using the volume element method as shown by Caillabet, et al. [START_REF] Caillabet | Computation of large-scale parameters for dispersion in fissured porous medium using finite-volume method[END_REF]. It is worth mentioning that Wood et al. [START_REF] Wood | Volume averaging for determining the e↵ective dispersion tensor: Closure using periodic unit cells and comparison with ensemble averaging[END_REF] proved that the closure problem solution for dispersion in heterogeneous porous media is equivalent to the one resulting from the ensemble averaging method.

The solution of the closure problem for homogeneous porous media was further investigated by Amaral Souto and Moyne [START_REF] Amaral-Souto | Dispersion in two-dimensional periodic porous media. Part II. Dispersion tensor[END_REF] considering ordered and disoredered geometries for the solid phase in the unit cell as well as the flow orientation for particle Reynolds number values ranging in the laminar-inertial regime.

The pertinence of using these geometries for the closure problem solution was validated by comparison with experiments in a subsequent work by Didierjean, et al. [START_REF] Didierjean | Dispersion in periodic porous media. Experience versus theory for two-dimensional systems[END_REF]. Later on, Wood [START_REF] Wood | Inertial e↵ects in dispersion in porous media[END_REF] studied the role of inertial e↵ects over dispersion in homogeneous porous media, finding that longitudinal dispersion was not dramatically a↵ected by inertia, whereas transverse dispersion was enhanced by a factor of 40 with respect to creeping flow conditions. Recently, Aguilar-Madera et al. [START_REF] Aguilar-Madera | E↵ective mass di↵usion and dispersion in random porous media[END_REF] reported that the flow direction is the main cause of anisotropy of the dispersion tensor, especially for the transverse component of the tensor.

Application of the volume averaging method has not been restricted to passive dispersion under one-phase flow condition; Quintard and Whitaker [START_REF] Quintard | Convection, dispersion, and interfacial transport of contaminants: Homogeneous porous media[END_REF] studied active dispersion, which corresponds to the dissolution of trapped non-aqueous phase liquids (NAPL) in the water phase that saturates homogeneous porous media.

These authors solved the corresponding closure problems in simple unit cells and later on Ahmadi et al. [START_REF] Ahmadi | Calculation of the e↵ective properties describing active dispersion in porous media: from simple to complex unit cells[END_REF] carried out the computations in complex unit cells involving thousands of pores.

In all the above applications of the volume averaging method, the closure problem formulation follows essentially the same philosophy outlined by Carbonell and Whitaker [START_REF] Carbonell | Dispersion in pulsed systems II: Theoretical developments for passive dispersion in porous media[END_REF]. In this approach, the closure problem requires knowledge of the pore-scale velocity fields and their deviations. In addition, the total dispersion tensor has been defined as the sum of an e↵ective di↵usivity tensor and a hydrodynamic dispersion tensor, with both contributions being dependent on the flow rate. This is inconsistent with the definition of the e↵ective di↵usivity, which is an intrinsic e↵ective medium coe cient that only depends of the porous medium geometry. Within this context, the purpose of the present work is twofold: firstly, we reformulate the closure problems in a way that they are only written in terms of closure variables, hence providing a fully upscaled closed form, and secondly, a new formulation of the total dispersion tensor, involving an e↵ective di↵usivity tensor and two hydrodynamic dispersion tensors is presented. The paper is organized as follows. In section 2, we provide some generalities about the upscaling process using the method of volume averaging and directly show the structure of the averaged model and the ancillary closure problems for mass and momentum transport, the latter corresponding to the inertial regime. In Section 3, we reformulate the e↵ective medium coe cients and the closure problems to meet our goals. Symmetry of the total dispersion tensor is then investigated. Results of the closure problems, which were solved in unit cells of two-dimensional model porous media for several flow rates and orientations are reported in Section 4. Finally, the ensuing conclusions are presented in Section 5.

Averaging

In order to upscale the pore-scale governing equations by means of the volume averaging method, it is necessary to define an averaging domain V of norm V that contains both solid and fluid phases such as the one sketched in Fig. 1. In terms of this averaging domain, let us define the superficial and intrinsic averaging operators as [START_REF] Whitaker | The Method of Volume Averaging[END_REF] h i = 1

V Z V dV (5a) h i = 1 V Z V dV (5b)
with being a piece-wise smooth function defined in the -phase occupying the region V of norm V within V . These two averaging operators are related by

h i = " h i (5c)
where " ⌘ V /V is the volume fraction of the fluid phase contained in the averaging domain, i.e. the porosity in the present development.

The averaging process, as reported in [START_REF] Whitaker | The Method of Volume Averaging[END_REF], commences with the application of the superficial averaging operator to the pore-scale equations. While doing this, it is necessary to interchange spatial di↵erentiation and integration and this is performed by means of the spatial averaging theorem [START_REF] Howes | The spatial averaging theorem revisited[END_REF] hr i = rh i + 1

V Z A n dA ( 6
)
where n is the unit normal vector at A pointing towards the -phase. The introduction of an interfacial integral resulting from the application of the spatial averaging theorem allows the interfacial boundary conditions of the pore-scale problem to be substituted when appropriate. However, at this stage, the resulting averaged equations contain both microscopic ( ) and macroscopic (h i) quantities. To progress towards a model involving macroscale variables only, it is necessary to express pointwise quantities in terms of their intrinsic averages and spatial deviations as proposed in [START_REF] Gray | A derivation of the equations for multiphase transport[END_REF] = h i + ˜ [START_REF] Aris | On the dispersion of a solute in a fluid flowing through a tube[END_REF] To carry out the development up to this stage where macroscopic equations are yet unclosed, it is not necessary to impose any length-scale constraint, as explained by Wood and Valdés-Parada [START_REF] Wood | Volume averaging: Local and nonlocal closures using a Green's function approach[END_REF]. Thus, the average model contains practically the same amount of information as its pore-scale counterpart. The systematic use of a set of time and length-scale constraints and assumptions in the form of scaling postulates is the essence of the upscaling procedure and allows filtering out unnecessary information at the macroscopic level. A careful and detailed explanation of the imposition of these postulates is available in ref. [START_REF] Wood | Volume averaging: Local and nonlocal closures using a Green's function approach[END_REF]. For the purposes of this work, it su ces to summarize that the characteristic size of the averaging domain, r 0 , is usually taken to be much larger than the characteristic size of the pores (` ) and, simultaneously, it must be much smaller than the characteristic length associated to the macroscale (L), i.e. [START_REF] Whitaker | The Method of Volume Averaging[END_REF] ` ⌧ r 0 ⌧ L

In addition, in many transport processes, there is a disparity between the characteristic time scales associated to pore-scale transport (say t ⇤ ) and to macroscopic transport (say t ⇤ h i ) so that t ⇤ ⌧ t ⇤ h i . Under these circumstances, averaged properties can be treated as if they were uniform, both in space and time, within the integration domains, equal to their values at the centroid. This is supported by the fact that non-local terms resulting from Taylor expansions of the averaged quantities away from the centroid can be neglected as far as the distances over which their successive gradients experience significant variations are much larger than r 0 (see [START_REF] Whitaker | The Method of Volume Averaging[END_REF], section 1.4.3 for instance). A direct corollary of this is the following average constraint for the deviations fields

h ˜ i = 0 (9)
The final step of the volume averaging method is to derive expressions that relate ˜ with h i and/or its derivatives, i.e., to close the macroscale model. The closure process can be summarized by the three following main steps:

• derive the governing equations for the deviations by subtracting the unclosed averaged equations from the pore-scale equations;

• simplify the problem by imposing reasonable constraints and assumptions on the basis of orders of magnitude analyses (i.e., scaling postulates). In particular, when the space and time constraints mentioned above are satisfied, a quasi-stationary closure problem can be obtained;

• formally solve the closure problem in terms of average quantities in simplified but still representative domains that represent the essential pore-scale geometry such as a periodic unit cell. This last step can be performed by using Green's functions as detailed by Wood and Valdés-Parada [START_REF] Wood | Volume averaging: Local and nonlocal closures using a Green's function approach[END_REF].

In this way, the spatial deviations can be related to average quantities and their derivatives by means of a linear superposition in terms of closure variables that can be shown to be integrals of the associated Green's functions. It is worth mentioning that without the separation of time and length scale constraints indicated above supporting the assumption that average quantities can be regarded as constants within the integration domain, it would not be possible to express the spatial deviations in the form of the superposition mentioned above.

These are the essential elements of the application of the volume averaging method to the pore-scale equations. A detailed description of the derivation of the upscaled model for mass and momentum transport is available elsewhere [cf . [START_REF] Carbonell | Dispersion in pulsed systems II: Theoretical developments for passive dispersion in porous media[END_REF][START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF][START_REF] Whitaker | The Method of Volume Averaging[END_REF]. The resulting macroscopic mass transport model, for " not necessarily constant, can be recalled as

" @hc A i @t + r • " hv i hc A i = r • " D ⇤ • rhc A i (10)
In Eq. [START_REF] Plumb | Dispersion in heterogeneous porous media: 1. Local volume averaging and largescale averaging[END_REF], the total dispersion tensor D ⇤ is defined by

D ⇤ = D ⇤ ef f + D = D 0 B @I + 1 V Z A n f dA 1 C A | {z } D ⇤ ef f hṽ f i | {z } D (11)
Here D ⇤ ef f and D denote an e↵ective di↵usivity tensor and a hydrodynamic dispersion tensor, respectively, while f is a closure variable vector defined by

cA = f • rhc A i (12)
that solves the following boundary-value problem

ṽ + v • rf = D r 2 f , in V (13a) n • rf = n , at A (13b) f (r + l i ) = f (r), i = 1, 2, 3 (13c)
hf i = 0 (13d)
As explained above, in order to arrive at the form of the closure problem reported in Eqs. [START_REF] Caillabet | Computation of large-scale parameters for dispersion in fissured porous medium using finite-volume method[END_REF], one also has to satisfy the time-scale constraint expressed by t t 0 with

1 ⌧ D t 0 `2 (14)
as reported in [START_REF] Whitaker | The Method of Volume Averaging[END_REF], section 3.3.4.

For momentum transport, the resulting equation from the upscaling process is the Darcy-Forchheimer equation [START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF]

hv i = K µ • rhp i F • hv i (15)
with K and F being the intrinsic permeability tensor and the Forchheimer correction tensor, respectively.

Note that this expression does not include the Brinkman correction term and it is thus constrained to the homogeneous portions of the system where this term can be shown to be negligible. Certainly, the above expression can be rearranged in a more compact form that resembles to Darcy's law:

hv i = H µ • rhp i (16)
where the tensor H may be regarded as an apparent permeability tensor, which is defined as

H 1 = K 1 • (I + F ) (17)
In order to compute this tensor, it is necessary to solve the following closure problem [START_REF] Lasseux | On the stationary macroscopic inertial e↵ects for one phase flow in ordered and disordered porous media[END_REF]:

r • E = 0, in V (18a) ✓ ⇢ v µ ◆ • rE = re + r 2 E + I, in V (18b) E = 0, at A (18c) E (r + l i ) = E (r), i = 1, 2, 3 (18d)
e (r + l i ) = e (r), i = 1, 2, 3 (18e)
he i = 0 (18f) hE i = H (18g)
where the vector e and the tensor E are the closure variables that are defined by [START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF] p = " µ e • (H )

1 • hv i = e • rhp i (19a) ṽ = ⇣ " E • (H ) 1 I ⌘ • hv i = 1 µ E + hE i • rhp i = 1 µ Ẽ • rhp i (19b)
with the obvious notation Ẽ = E hE i .

As mentioned above, the domain in which the closure problems are solved is usually a periodic unit cell (of period l i in the direction i = 1, 2, 3) that captures the essential features of the microscopic geometry, flow and mass transport, such as the one sketched in Figure 2 [START_REF] Amaral-Souto | Dispersion in two-dimensional periodic porous media. Part II. Dispersion tensor[END_REF].

Directing the attention to the closure problem for mass transport, it is worth noting that f depends on D , the fluid velocity within the pores, v , its deviations, and of the unit cell geometry. Similarly, the closure problem for momentum in Eqs. [START_REF] Aguilar-Madera | E↵ective mass di↵usion and dispersion in random porous media[END_REF] depends on the cell geometry, v and the fluid properties. Once the closure problem for f is solved, the fields are substituted into the integral of Eq. ( 11) to compute the total dispersion tensor. These integrals play the role of filters of redundant information coming from the closure problem (see [START_REF] Whitaker | The Method of Volume Averaging[END_REF]). In this way, the surface integral term in Eq. ( 11) is regarded as a di↵usive filter and it is contained in the e↵ective di↵usivity definition, whereas the volumetric integral in the third term is a convective filter and it constitutes the hydrodynamic dispersion term.

The structure of Eq. ( 11) nicely indicates that total dispersion is the sum of the e↵ective di↵usivity and hydrodynamic dispersion. However, as noted recently [26, see Section 6 therein], the di↵usive filter in Eq. ( 11) allows passing some convective information, (f is, in general, a function of the fluid velocity), so that

D ⇤
ef f changes with the fluid velocity, which is unexpected for a di↵usion coe cient. Furthermore, the fact that the closure problems' solutions require knowledge of the pointwise velocity field, seems inconsistent with the philosophy of the volume averaging method, in which the computation of e↵ective medium coe cients requires solving closure problems and not problems for the pointwise physical variables. Indeed, under their current forms, boundary value problems in Eqs ( 13) and ( 18) appear to remain unclosed. From a practical point of view, this means that the determination of D ⇤ requires a prior computation of the pore scale flow problem over the unit cell (see Section 2 in [START_REF] Carbonell | Dispersion in pulsed systems II: Theoretical developments for passive dispersion in porous media[END_REF]) for any combination of flow and/or mass transport regimes.

Therefore, an alternative formulation that does not involve the pore scale velocity is highly desirable. In the following, these issues are carefully addressed so that the e↵ective di↵usivity is defined as an intrinsic coe cient and the closure problems are written only in terms of closure variables.

E↵ective coe cients and closure problems

In this section, a new formulation of the dispersion tensor in terms of the intrinsic e↵ective di↵usivity tensor and the hydrodynamic dispersion is presented. Furthermore, in section 3.2, the associated closure problems for mass and momentum transport, including both viscous and inertial e↵ects, are reformulated using closure variables only. A symmetry analysis on D ⇤ is provided in section 3.3. Finally, simplifications in the case of creeping flow are presented in section 3.4.

Reformulation of the dispersion tensor

Let us start our derivations by noting that, in the absence of convection, Eqs. ( 13) reduce to

r 2 f 0 = 0, in V (20a) n • rf 0 = n , at A (20b) f 0 (r + l i ) = f 0 (r), i = 1, 2, 3 (20c)
hf 0 i = 0 (20d)
Here,

f 0 = f | v =0 is intrinsic, i.e.
is only a function of the unit cell geometry since it does not depend on D nor v . This closure problem corresponds to the one given by eqs. (1.4-58) in [START_REF] Whitaker | The Method of Volume Averaging[END_REF], the solution of which allows computing the intrinsic e↵ective di↵usivity tensor, D ef f given by [see Eq. (1.4-62) in 21]

D ef f = D 0 B @I + 1 V Z A n f 0 dA 1 C A (21)
In this way, after a simple algebraic manipulation, we may write Eq. ( 11) as follows

D ⇤ = D ef f + D h = D ef f + D V Z A n (f f 0 ) dA hṽ f i (22)
The last two terms in the above expression define the true hydrodynamic dispersion tensor, D h . Certainly, for highly convective conditions, i.e., for conditions in which

D V R A n (f f 0 ) dA ⌧ hf ṽ i , one expects D h ' D ' D ⇤ . Alternatively, when convective e↵ects are insignificant, D ⇤ ' D ⇤ ef f ' D ef f D h ' D .
However, for cases in which convective mass transfer is comparable to or greater than di↵usive transfer, one may expect di↵erent predictions for D h and D (as well as D ⇤ ef f being significantly di↵erent from D ef f ). This has been discussed in [26, see section 6]. As will be shown below, under these conditions, transverse dispersion may be larger than longitudinal dispersion.

At this point, it may appear that once the flow problem is solved over the unit cell, two closure problems need to be solved in order to compute D ⇤ . Actually, it is only necessary to solve the closure problem given by Eqs. [START_REF] Caillabet | Computation of large-scale parameters for dispersion in fissured porous medium using finite-volume method[END_REF], firstly under non-convective conditions and secondly for the desired flow condition. Let us now rewrite the closure problems and e↵ective medium coe cients for both momentum and mass transport in such a way that they do not require prior solution of the pore-scale flow problem.

Dispersion tensor determination from closure problems

We start with the reformulation of the closure problem for momentum given in Eqs. [START_REF] Aguilar-Madera | E↵ective mass di↵usion and dispersion in random porous media[END_REF]. Using the expressions for ṽ , hv i and H given in Eqs. (19b), ( 16) and (18g), yields v = 1 µ E • rhp i . As a consequence, the momentum-like equation (18b) can be rewritten as

⇢ µ 2 rhp i • E T • rE = re + r 2 E + I (23)
It is worth noting that, in the classical formulation, the pore-scale velocity must be solved every time the macroscopic pressure gradient changes, and in the current formulation, the link between v and rhp i is the closure variable, E which is finally integrated in the closure problem. In this way, one needs to specify 190 the macroscopic pressure gradient applied on the unit cell instead of providing the microscale velocity field.

Turning our attention to the mass transport closure problem, Eq. (13a) can be similarly reformulated as

µ 1 ⇣ Ẽ • rhp i + E • rhp i • rf ⌘ = D r 2 f (24)
In the same way, the hydrodynamic dispersion tensor takes the form

D h = D V Z A n (f f 0 ) dA + µ 1 rhp i • hE T f i (25)
While deriving this last expression, we have made use of the fact that, due to the constraint in [START_REF] Carbonell | Dispersion in pulsed systems II: Theoretical developments for passive dispersion in porous media[END_REF], rhp i can be assumed constant within the averaging domain and we have also considered Eq. (13d) which

195 implies h ẼT f i = hE T f i .
A convenient dimensionless form, well adapted for numerical solution, may now be used based on the following dimensionless quantities:

V ⇤ = V /`3, A ⇤ = A /`2, f ⇤ = f /`, f ⇤ 0 = f 0 /`, e ⇤ = e /`and E ⇤ = E /`2,
where `is the characteristic size of the unit cell (see figure 2). In addition, we introduce = rhp i /krhp i k the unit vector in the direction of the average pressure gradient. Using the same symbol r for the nabla operator with or without dimension (the latter being the product of the former by `) yields the following dimensionless closure problem for momentum transport

r • E ⇤ = 0, in V (26a) Re ⇤ • E ⇤T • rE ⇤ = re ⇤ + r 2 E ⇤ + I, in V (26b)
E ⇤ = 0, at A (26c) E ⇤ (r + l i ) = E ⇤ (r), i = 1, 2, 3 (26d)
e ⇤ (r

+ l i ) = e ⇤ (r), i = 1, 2, 3 (26e)
he ⇤ i = 0 (26f)
from which the dimensionless apparent permeability can be deduced as

H ⇤ = H /`2 = hE ⇤ i. In equation (26b
), Re ⇤ is the Reynolds number defined as Re ⇤ = ⇢ `3krhp i k/µ 2 . It must be noted that and Re ⇤ are the two parameters allowing to specify the pressure gradient direction and intensity, respectively.

Accordingly, the dimensionless closure problem for mass transport may be expressed as

P e khE ⇤ i • k • ⇣ Ẽ⇤T + E ⇤T • rf ⇤ ⌘ = r 2 f ⇤ (27a) n • rf ⇤ = n , at A (27b) f ⇤ (r + l i ) = f ⇤ (r), i = 1, 2, 3 (27c)
hf ⇤ i = 0 (27d)
Note that in Eq. (27a), the fields of E are now required in order to solve the closure problem for mass transport. They are obtained by solving Eqs. [START_REF] Valdés-Parada | The role of tortuosity in upscaling[END_REF], which, contrary to the original formulation [START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF], do not require the solution of the pore scale velocity and can now be solved as an independent boundary-value problem. In Eq. (27a), P e is the cell Péclet number defined as

P e = " khv i kD (28)
In this way, the total dispersion tensor can be expressed as

D ⇤ = D 0 B @I + 1 V ⇤ Z A n f ⇤ 0 dA ⇤ 1 C A | {z } e↵ective di↵usivity D ef f + D 0 B @ 1 V ⇤ Z A n f ⇤ f ⇤ 0 dA ⇤ + P e khE ⇤ i • k • hE ⇤T f ⇤ i 1 C A | {z } hydrodynamic dispersion D h (29)
where f ⇤ 0 is the solution of Eqs. [START_REF] Papathanasiou | A computational evaluation of the Ergun and Forchheimer equations for fibrous porous media[END_REF] when P e = 0. At this point of the analysis, the following comments are in order:

• The e↵ective di↵usion coe cient, D ef f , is only a function of the porous medium microstructure and it is no longer a↵ected by convective e↵ects. Although the diagonal terms of this tensor remain always positive, the surface integral part is a negative definite function that reflects the influence of the pore scale geometry.

• The components of the hydrodynamic dispersion tensor, D

h , only depend on: 1) The porous medium microstructure, 2) The cell Péclet number, P e, 3) The Reynolds number, Re ⇤ and 4) The macroscopic pressure gradient orientation defined by .

• Closure problems, as expressed in Eqs. ( 26) and [START_REF] Papathanasiou | A computational evaluation of the Ergun and Forchheimer equations for fibrous porous media[END_REF], are certainly di↵erent from their original formulations since they no longer require solving the pore scale flow field. In fact, for a given Re ⇤ and , the closure problem for total mass and momentum transport can be solved once in order to obtain the fields of E ⇤ so that D h can then be computed for any required value of the Péclet number after solving the closure problem for mass transport.

The forms of the closure problems given by Eqs. ( 26) and ( 27) are convenient if one is willing to compute D ⇤ for a prescribed macroscopic pressure gradient applied to the structure, as reflected in the expression of Eq. [START_REF] Auriault | On the asymmetry of the dispersion tensor in porous media[END_REF]. In some circumstances, it might also be of interest to determine D ⇤ for a specific mean flow direction that does not necessarily coincide with the corresponding pressure gradient. To do so, it is convenient to rewrite Eqs. (26b), (27a) and ( 29) in terms of the unit vector v in the direction of hv i, which is defined as

v = hv i khv ik = H ⇤ • kH ⇤ • k ( 30
)
When this is done, the momentum-like equation in the closure problem [START_REF] Valdés-Parada | The role of tortuosity in upscaling[END_REF] takes the form

Re ⇤ v • hE ⇤T i 1 khE ⇤ i 1 • v k • E ⇤T • rE ⇤ = re ⇤ + r 2 E ⇤ + I (31)
while the transport equation in the closure problem ( 27) is given by

P e v • hE ⇤T i 1 • ⇣ Ẽ⇤T + E ⇤T • rf ⇤ ⌘ = r 2 f ⇤ (32)
This yields the following expression for

D ⇤ D ⇤ D = I + 1 V ⇤ Z A n f ⇤ 0 dA ⇤ + 1 V ⇤ Z A n f ⇤ f ⇤ 0 dA ⇤ P e v • hE ⇤T i 1 • hE ⇤T f ⇤ i (33)
Equations ( 31) and (32) hold a stronger degree of non linearity than the versions written in terms of and certainly feature closure problems that are more complex to solve, although in Eq. (31), the prefactor 1/khE ⇤ i 1 • v k might be lumped together with Re ⇤ to form a cell Reynolds number given by

Re = Re ⇤ /khE ⇤ i 1 • v k = Re ⇤ khE ⇤ i • k = ⇢ khv ikμ .
While Re ⇤ and Re conveniently measure the inertial to viscous forces ratio at the scale, `, of the unit cell, for a prescribed pressure gradient and flow rate respectively, other definitions may be used when a di↵erent scaling is employed. In fact, one might use a Reynolds number based on p k, k being the magnitude of the intrinsic permeability tensor for an isotropic medium [START_REF] Lasseux | On the stationary macroscopic inertial e↵ects for one phase flow in ordered and disordered porous media[END_REF][START_REF] Papathanasiou | A computational evaluation of the Ergun and Forchheimer equations for fibrous porous media[END_REF], or, more commonly, a pore-scale Reynolds number, Re p =

⇢ khv i k` µ = Re ⇤ ` khE ⇤ i• k " `.
Accordingly, a pore Péclet number P e

p = khv i k` D = P e ` " ìs
often employed [START_REF] Amaral-Souto | Dispersion in two-dimensional periodic porous media. Part II. Dispersion tensor[END_REF][START_REF] Salles | Taylor dispersion in porous media. determination of the dispersion tensor[END_REF] and, for the sake of consistency with data reported in the literature, our results in section 4 will be presented in terms of Re p and P e p .

Symmetry properties of the dispersion tensor

Our purpose in this section is to investigate the symmetry properties of D ⇤ . The starting point of the development is to form the dyadic product of Eq. (27a) with f ⇤ which gives

P e khE ⇤ i • k • ⇣ hE ⇤T i f ⇤ E ⇤T • rf ⇤ + I f ⇤ ⌘ = r 2 f ⇤ f ⇤ (34)
The right hand side (rhs) of Eq. ( 34) can be rewritten in an equivalent form as

r 2 f ⇤ f ⇤ = r • rf ⇤ f ⇤ rf ⇤ T • rf ⇤ (35)
and while taking the dimensionless superficial average of Eq. (34) in which the rhs is replaced by its expression in Eq. ( 35), we obtain

P e khE ⇤ i • k • ⇣ hE ⇤T i hf ⇤ i hE ⇤T • rf ⇤ + I f ⇤ i ⌘ =r • h rf ⇤ f ⇤ i + 1 V ⇤ Z A n • rf ⇤ f ⇤ dA ⇤ h rf ⇤ T • rf ⇤ i (36)
Due to periodicity, the first term in the rhs of this last expression is such that r • h

⇣ rf ⇤ ⌘ f ⇤ i = 0 while the
second term in the rhs can be re-written as 1

V ⇤ R A n • ⇣ rf ⇤ ⌘ f ⇤ dA ⇤ = 1 V ⇤ R A n f ⇤
dA ⇤ upon making use of the boundary condition of Eq. (27b). Moreover, because of the zero average constraint of f ⇤ (see Eq (27d)), the first term in the left hand side of Eq. ( 36) is also zero, so that this equation can be written as 1

V ⇤ Z A n f ⇤ dA ⇤ P e khE ⇤ i • k • hE ⇤T f ⇤ i = P e khE ⇤ i • k • hE ⇤T • rf ⇤ f ⇤ i h rf ⇤ T • rf ⇤ i (37)
When this last result is introduced back into the expression of

D ⇤
D given by Eq. ( 29), we have

D ⇤ D = I h rf ⇤ T • rf ⇤ i + P e khE ⇤ i • k • hE ⇤T • rf ⇤ + 2I f ⇤ i (38)
While the first two tensors in the rhs of Eq. (38) are symmetric, the last term is not symmetric, except when the unit cell is symmetric and is along one of its symmetry axis. This features a non-symmetric D ⇤ in general. Moreover, by taking P e = 0, it straightforwardly follows from Eq. (38) that D ef f is a symmetric tensor. As a consequence, it is clear that asymmetry may only occur in D h . More specifically, a comparison between the expressions of D ⇤ in Eqs. ( 29)cand (38), taking Eq. (37) into account, indicates that asymmetry may originate from the volume average and/or the area integral parts of D h in the rhs of Eqs. [START_REF] Auriault | On the asymmetry of the dispersion tensor in porous media[END_REF]. A quick order of magnitude analysis to the formal solution of the closure problem given by eqs. [START_REF] Caillabet | Computation of large-scale parameters for dispersion in fissured porous medium using finite-volume method[END_REF] indicates that f = O (` ) and hence that the second term in the rhs of Eq. ( 38) is O (1). Consequently, a su cient condition for D ⇤ to be quasi-symmetric is when

O ✓ P e khE ⇤ i • k • hE ⇤T • rf ⇤ + 2I f ⇤ i ◆ ⌧ 1 (39)
Since rf ⇤ can be estimated to be O (1), the above constraint for D ⇤ to be quasi-symmetric is finally

P e ⌧ O ✓ " ` ◆ (40)
or, equivalently, in terms of the pore Péclet number, P e p = khv i k` D P e p ⌧ 1 (41) in agreement with some investigations reported in [START_REF] Auriault | On the asymmetry of the dispersion tensor in porous media[END_REF]. The su cient constraint given above is independent of the Reynolds number, and this is an important feature of the symmetry properties of the dispersion tensor that requires a closer attention regarding the compatibility of symmetry with the existence of inertial e↵ects.

Since inertia becomes significant when Re p & 10, the su cient condition (41) for symmetry would require, in that case, that the Schmidt number, Sc, is such that

Sc = P e p Re p = µ ⇢ D ⌧ 0.1 ( 42
)
This condition is never met for conventional fluids as Sc = O(1) for gases while, for liquids, it is rather O( 103 ). For porous structures in which the unit cell does not possess any specific geometrical symmetry or when the mean flow is not along a symmetry axis, this suggests that inertia is a potential mechanism that can certainly trigger asymmetry of D ⇤ .

Finally, one should note that the condition given in (41) is independent from (but compatible with) the time-scale constraint expressed in ( 14) that is required to treat the mass transport closure problem as quasi-steady.

Creeping flow regime

Special attention should be dedicated to situations in which the Reynolds number is such that the flow remains in the creeping regime, i.e. typically when Re p ⌧ 10. In such circumstances, the closure problem for momentum transport does not depend on the Reynolds number nor on the macroscopic pressure gradient orientation and is hence intrinsic, yielding H = K . As a consequence, the computation of D ⇤ requires the momentum closure problem to be solved only once and the solution can be used for any Péclet number value and .

In the following sections, some computational results are provided in both the inertial and creeping flow regimes.

Results

The dimensionless closure problems given by Eqs. ( 26) and ( 27) were solved on periodic two-dimensional unit cells of model porous media using the finite element software Comsol Multiphysics 4.4 involving su cient mesh elements to guarantee consistency in the numerical results. A mesh-convergence analysis showed that a dimensionless grid-block area of about 4 ⇥ 10 5 along with a dimensionless boundary element size of 0.01 were appropriate to achieve convergent numerical results. The solid phase was considered as parallel cylinders of square cross sections arranged on a square regular periodic pattern (see figure 2). Two di↵erent values of the porosity, namely " = 0.8 and " = 0.4, were considered. Dispersion is studied in the plane orthogonal to cylinders axes as a function of P e p , without inertia (Re p = 0) and for Re p = 200 while considering ✓ = 0, ⇡/8 and ⇡/4. For the pore-length ` in the definition of Re p and P e p a quantitative estimate, based on the hydraulic diameter and already proposed in [START_REF] Eidsath | Dispersion in pulsed systems III: Comparison between theory and experiments for packed beds[END_REF] (see also [START_REF] Whitaker | The Method of Volume Averaging[END_REF] p. 144), is employed as ` ⇠ ` " 1 " , ` being the cylinder size (see figure 2). This yields

Re p = Re ⇤ khE ⇤ i• k` (1 " )`a nd P e p = P e ` ( 1
" )2
Although closure problems are solved for a prescribed pressure gradient direction defined by , inclined of an angle ✓ on e x , it is more physically appealing to present the results in terms of the orientation of the average velocity, i.e. in the longitudinal and transverse directions v (see Eq. ( 30)) and ⌧ v , ⌧ v being the unit vector directly orthogonal to v . Note that, even for porous materials such that K is a spherical tensor, as for the structure under consideration here, and v are not necessarily aligned when Re p 6 = 0

(see [START_REF] Lasseux | On the stationary macroscopic inertial e↵ects for one phase flow in ordered and disordered porous media[END_REF]) and we denote by ✓ v the macroscopic flow inclination on the horizontal x-direction. As a matter of fact, ✓ = ✓ v for ✓ = 0 and ⇡/4; however, for ✓ = ⇡/8 we obtained that ✓ v = ⇡/7 for " = 0.4 and ✓ v ⇡ 0.1521 for " = 0.8, in both cases Re p = 200. The longitudinal (along v ) and transverse (along ⌧ v ) components of the total dispersion tensor can be obtained from the components of D ⇤ in (e

x , e y ) as follows,

D ⇤ , v v = cos 2 ✓ v D ⇤ ,xx + sin 2 ✓ v D ⇤ ,yy + sin ✓ v cos ✓ v (D ⇤ ,yx + D ⇤ ,xy ) (43a) D ⇤ ,⌧v⌧v = sin 2 ✓ v D ⇤ ,xx + cos 2 ✓ v D ⇤ ,yy sin ✓ v cos ✓ v (D ⇤ ,yx + D ⇤ ,xy ) (43b)
When mass transport is strongly convection-dominated (i.e. when the Péclet number is large compared to unity), D ⇤ essentially depends on the volume average part of the hydrodynamic dispersion represented by the last term in Eq. ( 29) (or (33)). Examples of the fields of the v v component of the corresponding tensor

•E ⇤T f ⇤ khE ⇤ i• k
are reported over a unit cell in figures 3 and 4 for P e p = 10 3 , the three flow angles and the two Reynolds numbers Re p = 0 and 200. Clearly, the flow orientation and the inertial e↵ects have a determinant e↵ect on these fields and consequently on the components of the total dispersion tensor as will be shown below. Although strongly correlated to the flow structure as indicated by the flow streamlines that are superimposed in Figs. 3 and4, the dependence of the fields on Re p and ✓ is complex as their magnitude varies non-monotonically with ✓ and increases or decreases with Re p . This is reflected on the graphs of Fig. 5 and6 where the dimensionless longitudinal and transverse dispersion coe cients are represented versus

P e

p for the di↵erent values of ✓ and Re p . Regarding these results, the following comments are in order: 1. For P e p  1, all the predictions of D ⇤ converge to the same value that is given by the e↵ective di↵usivity, which is insensitive to inertial e↵ects and to the flow angle. As a matter of fact, under these conditions, D ⇤ reduces to D ef f , which only requires solving the closure problem for f 0 , given by eqs. (20).

For P e

p > 1, the longitudinal dispersion increases when Re p increases, whatever the flow angle, for a su ciently large porosity (see Figs. 5i-iii for " = 0.8). However, for a smaller porosity (" = 0.4, see Figs. 5a-c), the longitudinal dispersion may decrease or increase with Re p . In particular, for P e p & 100, D ⇤ , v v increases in the presence of inertia when the flow is not aligned with the principal axes (see Fig. 5b andc) but decreases with Re p when ✓ = 0 (Fig. 5a). This surprising behavior contrasts with the case " = 0.8 for the same range of P e p and the same flow angle (see Fig. 5i). A physical explanation can be given from the fields in Figs. 3a, 3i and4a, 4i). On the one hand, for ✓ = 0, Re p = 0 and " = 0.4, the vertical gaps between two adjacent cylinders in the direction of the flow are occupied by vortices taking place at small velocities. On the other hand, for " = 0.8 and the same flow conditions, only two small vortices are present in these regions that are mainly filled with tortuous flow streamlines connected from the entrance to the exit of the unit cell. When Re p = 200 and the same flow angle, the vortices pattern is not significantly modified in the unit cell with " = 0.4 (only two additional vortices are produced and velocities remain extremely small) whereas streamlines are significantly straightened in the horizontal channels in the vicinity of the vertical gaps. The conjunction of these two coupled mechanisms yields a less e cient dispersion and a decrease of D ⇤ , v v in the presence of inertia for ✓ = 0 and " = 0.4. Conversely, for " = 0.8, when Re p = 200 (✓ = 0), large eddies with significant velocities are taking place in the vertical gaps. Although they induce straighter streamlines in the horizontal channels where spreading is however favored, the net result is a significant improvement of dispersion and hence a larger value of D ⇤ , v v . These last two mechanisms, i.e. development of eddies with significant velocities allowing the development of spreading, are the main ones explaining the increase of the longitudinal dispersion with Re p for all other flow angles and the two values of " . 3. Keeping all other parameters the same, the longitudinal dispersion coe cient increases when porosity decreases. Moreover, the exponent of the power-law dependence of D ⇤ , v v on P e p , occurring at su ciently large values of the Péclet number, is not modified by the presence of inertia, which however favors its emergence at smaller P e p . 4. The transverse dispersion coe cient also exhibits a complex dependence upon P e p , ✓ and Re p . As for longitudinal dispersion, the transverse coe cient increases when porosity decreases for a given P e p , ✓ and Re p .

In figures 7, we have reported the ratio D ⇤ , v v /D ⇤ ,⌧v⌧v versus P e p for the values of " , ✓ and Re p under consideration. As expected, for P e p  1, mass transport is mainly driven by di↵usion and consequently

D ⇤ , v v = D ⇤ ,⌧v⌧v
. Interestingly, for 1 . P e p . 100 and ✓ = ⇡/4, the transverse component is larger than its longitudinal counterpart, except for " = 0.8, Re p = 200. This observation is consistent with numerical results from Salles et al. [START_REF] Salles | Taylor dispersion in porous media. determination of the dispersion tensor[END_REF] (see Table 5 therein). However, this e↵ect is no longer present when flow is along the principal axes of the structure and for a large value of the porosity when inertia is significant. Certainly, for P e p larger than ⇠ 100, the longitudinal component is larger than the transverse component, regardless of ✓ and the flow regime. Under these conditions, the longitudinal to transverse dispersion coe cients ratio corresponding to ✓ = ⇡/4 is larger than the values corresponding to ✓ = ⇡/8 for cases involving (or not) inertial e↵ects and for the two values of " . These observations are in agreement with a previous study [see,

Conclusions

In this work, the e↵ective di↵usivity and hydrodynamic dispersion tensors that compose the total dispersion tensor have been redefined in a volume averaging context. In these reformulations, the e↵ective di↵usivity is no longer a function of the flow and the ancillary closure problem corresponds to the one typically found when studying passive di↵usion in porous media. The hydrodynamic dispersion tensor was shown to depend on the Péclet number, the Reynolds number, the macroscopic flow angle and the geometry of the porous medium represented by the unit cell. It was shown that the dispersion tensor can be obtained from the solution of two closure problems, namely, one for momentum and one for mass transport. Both closure problems were formulated in terms of closure variables only, showing that the pore scale flow fields are not necessary. In this formulation, the mass transport closure problem needs to be solved twice, firstly with a zero Péclet number yielding the e↵ective di↵usivity and secondly for the desired Péclet number and flow angle. In the latter case, the solution of the closure problem for momentum, obtained with the required Reynolds number and flow angle, is needed. In the creeping flow regime, the closure problem for momentum is intrinsic and needs to be solved only once. This is a clear di↵erence compared to the existing formulation in which the pore scale flow problem was a necessary input for both the momentum and mass transport closure problems. The benefits of this new formulation are at least twofold. Firstly, it provides a completely closed set of boundary value problems as intended in the volume averaging method. Secondly, since the closure problem for mass transport is now coupled to the one for momentum transport, for conditions in which Re p 1, the numerical solution of the latter is significantly less demanding in terms of computational resources in the current formulation. Symmetry of the total dispersion tensor was analyzed and it was shown that D ⇤ is not symmetric in the general case. The condition for D ⇤ to be quasi-symmetric was investigated and a su cient condition is when the pore Péclet number remains small compared to unity, a situation that is however not compatible with the existence of significant inertial e↵ects for conventional fluids.

The ancillary closure problems were solved on periodic unit cells of 2D model structures made of parallel cylinders of square cross section arranged on a square regular pattern. The closure problem solutions were used to compute the longitudinal and transverse components of the total dispersion tensor. This was carried out considering three flow angles, with and without inertial e↵ects. For the mass transport and flow conditions under consideration, the dispersion coe cients are highly sensitive to the flow rate and to the flow angle. Furthermore, it was shown that, for small enough values of the porosity, inertial e↵ects may lead to a smaller longitudinal dispersion coe cient compared to that in creeping flow conditions, specifically when flow is along a principal axis of the structure. For some particular flow orientations within the structure under consideration (namely ✓ = ⇡/4), and except for a large value of the porosity when inertia is significant, the transverse component may be larger than the longitudinal dispersion for 1 . P e p . 100. For P e p & 100, the longitudinal component is larger than the transverse component for any flow condition and porosity.

In essence, the methodology developed in the present work could be extended to broader situations such as heat transfer in homogeneous porous media involving convection in the fluid phase as studied in [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: Numerical results for e↵ective properties[END_REF] or dispersion in heterogeneous porous media investigated in [START_REF] Quintard | Dispersion in heterogeneous porous media: One-equation nonequilibrium model[END_REF], among many others. In both cases, the transport equations that involve di↵usion and convective mechanisms are formally the same as Eq.

(1). Therefore, it is not hard to realize that the formulation of closure problems provided here can be straightforwardly applied to these and other similar situations.

E ⇤ (= E /`2) dimensionless version of E f closure variable that maps rhc A i onto cA , m f ⇤ (= f /`) dimensionless form of f f 0 (= f | v =0 ) closure

Figs. 2

 2 Figs. 2 and 3 in ref. 15].

A

  norm of the solid-fluid interface contained in the averaging domain, m 2 A ⇤ (= A /`2 ) dimensionless area of the solid-fluid interface c A molar concentration of species A in the -phase, mol/m 3 hc A i intrinsic average of c A , mol/m 3 cA spatial deviations of c A with respect to hc A i , mol/m 3 D species mixture di↵usivity in the -phase, m 2 /s D ⇤ total dispersion tensor, m 2 /s D hydrodynamic dispersion tensor, m 2 /s D ef f intrinsic e↵ective di↵usion tensor, m 2 /s D ⇤ ef f e↵ective di↵usion tensor that depends of the flow rate, m 2 /s e closure variable that maps rhp i onto p , m e ⇤ (= e /`) dimensionless version of e E closure variable associated to the spatial variations of the velocity, m 2

3 V 3 Figures captions Figure 1 .Figure 2 .

 33captions12 Figures captions

Figure 3 .Figure 4 .

 34 Figure 3. Examples of the fields of the v v component of the tensor • E ⇤T f ⇤ /kH • k and velocity streamlines for three flow angles taking " = 0.4 and P e p = 1000 for a) ✓ = 0, Re p = 0; b) ✓ = ⇡/8, Re p = 0; c) ✓ = ⇡/4, Re p = 0; i) ✓ = 0, Re p = 200; ii) ✓ = ⇡/8, Re p = 200; iii) ✓ = ⇡/4, Re p = 200.

Figure 5 .

 5 Figure 5. Longitudinal component of the total dispersion tensor vs. P e p taking three di↵erent values of ✓. Results are obtained from solving the associated closure problems in a 2D unit cell with the solid

Figure 6 .

 6 Figure 6. Transverse component of the total dispersion tensor vs. P e p taking three di↵erent values of ✓. Results are obtained from solving the associated closure problems in a 2D unit cell with the solid phase modeled as a square having porosities of 0.4 and 0.8. Black curves correspond to Re p = 0 and blue curves correspond to Re p = 200.

Figure 7 .

 7 Figure 7. Dependence of the ratio of the longitudinal and transverse components of the total dispersion tensor with P e p for three flow orientations taking " = 0.4 and " = 0.8. The results were obtained from solving the closure problems in unit cells with the solid phase represented by a square obstacle.

Figure 1 :Figure 2 :

 12 Figure 1: Sketch of the system, including characteristic lengths and a sample of the averaging domain.

Figure 3 :

 3 Figure 3: Examples of the fields of the v v component of the tensor • E ⇤T f ⇤ /kH • k and velocity streamlines for three flow angles taking " = 0.4 and P e p = 1000 for a) ✓ = 0, Re p = 0; b) ✓ = ⇡/8, Re p = 0; c) ✓ = ⇡/4, Re p = 0; i) ✓ = 0, Re p = 200; ii) ✓ = ⇡/8, Re p = 200; iii) ✓ = ⇡/4, Re p = 200.

Figure 4 :

 4 Figure 4: Examples of the fields of the v v component of the tensor • E ⇤T f ⇤ /kH • k and velocity streamlines for three flow angles taking " = 0.8 and P e p = 1000 for a) ✓ = 0, Re p = 0; b) ✓ = ⇡/8, Re p = 0; c) ✓ = ⇡/4, Re p = 0; i) ✓ = 0, Re p = 200; ii) ✓ = ⇡/8, Re p = 200; iii) ✓ = ⇡/4, Re p = 200.

Figure 5 :

 5 Figure 5: Longitudinal component of the total dispersion tensor vs. P e p taking three di↵erent values of ✓. Results are obtained from solving the associated closure problems in a 2D unit cell having porosities of 0.4 and 0.8 with the solid phase modeled as a square. Black curves correspond to Re p = 0 and blue curves correspond to Re p = 200.

Figure 6 :Figure 7 :

 67 Figure 6: Transverse component of the total dispersion tensor vs. P e p taking three di↵erent values of ✓. Results are obtained from solving the associated closure problems in a 2D unit cell having porosities of 0.4 and 0.8 with the solid phase modeled as a square. Black curves correspond to Re p = 0 and blue curves correspond to Re p = 200.

Acknowledgments

FJVP is grateful to École CentraleSupélec, Châtenay-Malabry and to the Ecole Nationale Supérieure des Arts et Métiers for the visiting professor positions. The same author is thankful to Fondo Sectorial de Investigación para la Educación from CONACyT (Project number: 12511908; Arrangement number: 112087) for the financial aid provided.