D Lasseux

F J Vald´es-Parada

Symmetry properties of macroscopic transport coefficients in porous media

We report on symmetry properties of tensorial effective transport coefficients characteristic of many transport phenomena in porous systems at the macroscopic scale. The effective coefficients in the macroscopic models (derived by upscaling (volume averaging) the governing equations at the underlying scale) are obtained from the solution of closure problems that allow passing the information from the lower to the upper scale. The symmetry properties of the macroscopic coefficients are identified from a formal analysis of the closure problems and this is illustrated for several different physical mechanisms, namely, one-phase flow in homogeneous porous media involving inertial effects, slip flow in the creeping regime, momentum transport in a fracture relying on the Reynolds model including slip effects, single-phase flow in heterogeneous porous media embedding a porous matrix and a clear fluid region, two-phase momentum transport in homogeneous porous media, as well as dispersive heat and mass transport. The results from the analysis of these study cases are summarized as follows. For inertial single-phase flow, the apparent permeability tensor is irreducibly decomposed into its symmetric (viscous) and skew-symmetric (inertial) parts; for creeping slip-flow, the apparent permeability tensor is not symmetric; for one-phase slightly compressible gas flow in the slip regime within a fracture, the effective transmissivity tensor is symmetric, a result that remains valid in the absence of slip; for creeping one-phase flow in heterogeneous media, the permeability tensor is symmetric; for two-phase flow, we found the dominant permeability tensors to be symmetric, whereas the coupling tensors do not exhibit any special symmetry property; finally for dispersive heat transfer, the thermal conductivity tensors include a symmetric and a skew-symmetric part, the latter being a consequence of convective transport only. A similar result is achieved for mass dispersion. Beyond the physical mechanisms under consideration in the present work, the reported technique can be viewed as a general methodology applicable to any type of upscaled model obtained by volume averaging.

I. INTRODUCTION

Modeling transport phenomena in porous media at the pore-scale is often not sufficient for applications as a macroscopic description of the physical mechanisms through fluxto-force relationships is of more practical interest. Derivation of macroscale models from the underlying physics at the pore scale is hence of major importance, and this can be carried out with the aid of an upscaling technique such as volume averaging, [START_REF] Whitaker | The Method of Volume Averaging[END_REF] homogenization, [START_REF] Auriault | Homogenization of Coupled Phenomena in Heterogeneous Porous Media[END_REF] thermodynamically constrained averaging theory, [START_REF] Gray | Introduction to the Thermodynamically Constrained Averaging Theory for Porous Medium Systems[END_REF] or stochastic approaches. [START_REF] Gelhar | Stochastic Subsurface Hydrology[END_REF] For most of them, the idea consists in the derivation of an upscaled model that involves macroscopic coefficients together with the appropriate boundary-value problems allowing the prediction of these coefficients for a given microstructure and transport conditions taking place in the porous medium. To complete the description and the understanding of the physics at the macroscopic scale, knowledge of the symmetry properties of the macroscopic transport coefficients is very important. In a general a) Email: didier.lasseux@u-bordeaux.fr sense, a tensorial transport coefficient relates the flux to the force in the macroscopic model. When this coefficient is symmetric, it indicates that the flux in the ith-direction due to a force in the jth-direction is exactly the same as the flux in the jth-direction due to the same force in the ith-direction, whatever the microstructure of the medium, and this represents fundamental information on the transport process. In other words, accounting for the symmetry properties of the coefficients involved in macroscale models is essential from a mathematical viewpoint, especially in terms of the manipulations related to the macroscale flux terms. The importance of the symmetry properties is also evident for the characterization of transport phenomena in multiscale systems from an experimental or numerical point of view, as the knowledge of the number of independent components of the tensors that need to be determined is of prime interest. Moreover, if an effective-medium equation is taken as the departing model for a subsequent upscaling process, knowledge of the symmetry properties of the associated coefficients is crucial to analyze the properties of the coefficients involved in the resulting largerscale model. A clear example of the above is the study of momentum transport in porous media. Since the derivation DOI : 10.1063/1.4979907 of the vector form of Darcy's law involves a permeability tensor, the properties of this tensor appear to be ambiguous when this macroscale model is applied to situations different from steady and incompressible creeping flow in rigid and homogeneous porous media, where it is well-known that this tensor is symmetric as will be recalled below. Under inertial flow conditions, there may be situations in which this tensor is quasi-symmetric or even skew-symmetric, depending on the value of the Reynolds number. If more than one phase is involved, the porous medium structure is heterogeneous, or the flow is not incompressible, the properties of the permeability tensor become more obscure. Certainly, these doubts are extensible to the dispersion tensor involved in heat or mass transport in porous media.

To be more precise, a quick review of the literature dealing with the analysis of the symmetry of macroscopic transport coefficients shows that this property has been reported for few transport mechanisms only. The intrinsic permeability tensor for one-phase incompressible momentum transport in the creeping regime is probably the one for which symmetry has been studied most frequently, from energy considerations, [START_REF] Matheron | Elements pour une Théorie des Milieux Poreux[END_REF] double-scale homogenization, [START_REF] Sanchez-Palencia | Nonhomogeneous Media and Vibration Theory[END_REF][START_REF] Auriault | Dynamics of porous saturated media, checking of the generalized law of Darcy[END_REF][START_REF] Skjetne | Homogenization of wall-slip gas flow through porous media[END_REF][START_REF] Auriault | Transport in porous media: Upscaling by multiscale asymptotic expansions[END_REF] stochastic approach, [START_REF] Dagan | Flow and Transport in Porous Formations[END_REF] or volume averaging. [START_REF] Whitaker | The Method of Volume Averaging[END_REF][START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF][START_REF] Lasseux | An improved macroscale model for gas slip flow in porous media[END_REF] When inertia is present, there is no formal result available regarding the symmetry properties of the apparent permeability tensor except some observation that it is not symmetric in general as obtained from the volume averaging procedure applied to the incompressible Navier-Stokes equations [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF] or from numerical simulations. [START_REF] Koch | Moderate Reynolds number flows through periodic and random arrays of aligned cylinders[END_REF][START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF] Nonsymmetry of the apparent permeability tensor for single-phase slip-flow was only mentioned in the work by Skjetne and Auriault [START_REF] Skjetne | Homogenization of wall-slip gas flow through porous media[END_REF] using homogenization. It was recently addressed in more detail by Lasseux et al. [START_REF] Lasseux | An improved macroscale model for gas slip flow in porous media[END_REF] by means of the volume averaging method. In the case of two-phase flow in homogeneous porous media, symmetry of the effective (or dominant) permeability tensors was reported together with reciprocity relationships for the coupling permeability tensors in accordance with Onsager's principle, both in the context of homogenization [START_REF] Auriault | Transport in porous media: Upscaling by multiscale asymptotic expansions[END_REF][START_REF] Auriault | Nonsaturated deformable porous media: Quasistatics[END_REF] and volume averaging. [START_REF] Whitaker | The closure problem for two-phase flow in homogeneous porous media[END_REF][START_REF] Lasseux | Determination of permeability tensors for two-phase flow in homogeneous porous media: Theory[END_REF] The case of mass dispersion of a tracer is strikingly a situation where confusion on the symmetry properties of the dispersion tensor reported in the literature has lasted for several decades. The dispersion tensor was sometimes assumed (without any proof) to be non-symmetric in the general case. [START_REF] Whitaker | The Method of Volume Averaging[END_REF][START_REF] Carbonell | Dispersion in pulsed systemsII: Theoretical developments for passive dispersion in porous media[END_REF][START_REF] Auriault | Taylor dispersion in porous media: Analysis by multiple scale expansions[END_REF][START_REF] Amaral-Souto | Dispersion in two-dimensional periodic porous media. Part I. Hydrodynamics[END_REF] Some proof was provided while employing double-scale homogenization [START_REF] Auriault | On the asymmetry of the dispersion tensor in porous media[END_REF] or volume averaging. [START_REF] Valdés-Parada | A new formulation of the dispersion tensor in homogeneous porous media[END_REF] Conversely, symmetry of this tensor is often implicitly postulated, [START_REF] Fel | Dispersion and dispersivity tensors in saturated porous media with uniaxial symmetry[END_REF] demonstrated when determined from the method of moments, [START_REF] Brenner | Dispersion resulting from flow through spatially periodic porous media[END_REF][START_REF] Salles | Taylor dispersion in porous media. Determination of the dispersion tensor[END_REF][START_REF] Koch | The symmetry properties of the effective diffusivity tensor in anisotropic porous media[END_REF] or simply considered as such because any nonsymmetric component in this tensor would be unimportant in the estimation of the dispersive flux in the macroscopic dispersion equation. [START_REF] Whitaker | The Method of Volume Averaging[END_REF][START_REF] Carbonell | Dispersion in pulsed systemsII: Theoretical developments for passive dispersion in porous media[END_REF][START_REF] Amaral-Souto | Dispersion in two-dimensional periodic porous media. Part I. Hydrodynamics[END_REF] It is sometimes invoked on the basis of Onsager's reciprocity arguments. [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Bear | Effects of material symmetry on the coefficients of transport in anisotropic porous media[END_REF] Onsager's arguments have been widely employed in the upscaling process from the molecular scale to the continuum scale to justify reciprocity of the coupling at the continuum scale. [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF][START_REF] Casimir | On onsager's principle of microscopic reversibility[END_REF] However, they do not necessarily apply for subsequent upscaling procedures leading to macroscopic transport models in porous media. [START_REF] Auriault | On the asymmetry of the dispersion tensor in porous media[END_REF][START_REF] Auriault | On the cross-effects of coupled macroscopic transport equations in porous media[END_REF] In all cases mentioned above, it must be noted that the symmetry analysis of macroscopic tensorial transport coefficients does not rely on a well-defined and systematic methodology. There are many other transport processes, apart from those cited above, for which, to the best of our knowledge, no results at all were reported so far on the symmetry properties of the macroscopic coefficients present in the upscaled models. The case of single-phase incompressible flow through fractures featuring heterogeneous aperture fields, or through strongly heterogeneous porous materials embedding porous matrix and clear fluid regions, together with heat dispersion are amongst the principal ones.

The focus of the present work is hence laid upon an efficient method dedicated to the analysis of the symmetry properties of macroscopic coefficients arising in upscaled transport models within the framework of the volume averaging method (detailed in Appendix A). With this purpose in mind, a series of problems dealing with momentum, heat, and mass transport in porous media are analyzed. In all cases, we limit our presentation to the pore-scale model, the upscaled model, and the ancillary closure problems, so that attention can be directed to the reformulation and symmetry properties analysis of the corresponding effective-medium coefficients.

The paper is organized as follows. In Section II, study cases related to momentum transport are considered, namely, inertial one-phase flow, slip flow in the creeping regime, flow in a fracture including slip effects, one-phase flow in a medium embedding a porous matrix and a clear fluid region, and finally creeping two-phase flow. In Section III, dispersive heat and mass transport in homogeneous porous media are subsequently analyzed, highlighting the relevance of the symmetric and skew-symmetric parts of the associated macroscopic transport coefficients. Finally, the corresponding conclusions are presented in Section IV. The derivations in all the study cases require the use of tensor analysis; thus for the sake of brevity in presentation, a list of the most relevant tensor identities is provided in Appendix B, together with some important recalls on vector and tensor algebra necessary to recover them.

II. MOMENTUM TRANSPORT

A. Inertial one-phase flow

Consider the steady, incompressible, and Newtonian flow of a single fluid phase β saturating a rigid and homogeneous porous medium, for which the solid matrix is referred to as the σ-phase. The configuration is schematized in Fig. 1 where we have represented the macroscopic region V M of scale L and the averaging domain V of radius r 0 over which the pore-scale physical model is averaged and where the βphase and σ-phase can be identified. The governing equations for total mass and momentum transport at the pore-scale are given by the density and dynamic viscosity of the β-phase for single phase-flow. As detailed by Whitaker, [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF] the resulting upscaled model arising from the volume averaging method consists of the macroscale continuity equation,

∇ • v β = 0, in V β , ( 1a
)
ρv β • ∇v β = -∇p β + µ∇ 2 v β , in V β , (1b) v β = 0, at A βσ . ( 1c
∇ • v β β = 0 (2a)
and the following Darcy-like equation:

v β = - H µ • ∇ p β β . ( 2b
)
Here the superficial averaged velocity and pressure are defined as (see Appendix A)

v β = 1 V V β v β dV , p β β = 1 V β V β p β dV . (3)
Note that v β = ε v β β , with ε being the porous medium porosity. In addition, the apparent permeability tensor H is defined as

H = E . ( 4
)
In the above expression, the second-order tensor E is a closure variable that maps v β β onto the spatial variations of the fluid velocity, ṽβ , and it satisfies the following boundaryvalue problem (see the work of Whitaker 11 for details in the derivation):

∇ • E = 0, in V β , (5a) ρ µ v β • ∇E = -∇e + ∇ 2 E + I, in V β , (5b)
E = 0, at A βσ , (5c) ψ(r) = ψ(r + l i ), ψ = E, e, (5d)
e β = 0. (5e)
At this point, it is worth mentioning that the developments presented this far involve adopting a set of length-scale constraints and assumptions as explained in item 5 of Appendix A. One of these scaling postulates is that a periodic unit cell is a reasonable solution domain for the above closure problem and thus the periodic boundary condition given in Eq. (5d). In this problem, the vector e is another closure variable that maps v β β onto the spatial variations of the fluid pressure, pβ . Note that ∇E is a third-order tensor which is defined in Eq. (B3) of Appendix B. This closure problem has been numerically solved by Lasseux et al. [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF] on unit cells of model structures for several Reynolds number values and incident flow angle orientations.

Our objective in this section is to analyze the structure of the apparent permeability tensor H. From Eq. ( 4), it follows that, in order to gain more knowledge about this coefficient, it is necessary to have more insight into the closure variable E. With this in mind, let us pre-multiply Eq. ( 5b) by E T and take into account the solenoidal nature of the fluid velocity fields, in order to obtain

ρ µ E T • ∇ • v β E = -E T • ∇e + E T • ∇ 2 E + E T . ( 6
)
Applying the superficial average operator to this equation leads to

ρ µ E T • ∇ • v β E = -E T • ∇e + E T • ∇ 2 E + H T , ( 7
)
where, for the last term of the above equation, we used the definition given in Eq. ( 4). Moreover, employing the identity in Eq. (B18) of Appendix B, taking into account Eq. (5a), and after using the spatial averaging theorem in its divergence form as indicated in Eq. (A2) of Appendix A, the first term on the right-hand side (rhs) of the above equation takes the form

E T • ∇e = ∇ • Ee = ∇ • Ee + 1 V A β σ n • EedA, (8)
where n denotes the unit normal vector to the A βσ interface directed from the βto the σ-phase (see Fig. 1). Note also that Ee is a third-order tensor resulting from the outer product defined in Eq. (B11) of Appendix B. When the interfacial boundary condition given by Eq. (5c) is considered and while limiting the analysis to the bulk of the porous medium where average quantities are position-invariant, it results that E T • ∇e = 0 and Eq. ( 7) can be written as

ρ µ E T • ∇ • v β E = E T • ∇ 2 E + H T . ( 9
)
Directing the attention to the first term on the rhs of this equation, and using successively the identities reported in Appendix B in Eqs. (B21), (B12), and (B22) (taking A to be E T and B to be E) and Eq. (B16), we have

E T • ∇ 2 E = ∇ 2 E T • E T = ∇ • ∇E T • E T -(∇E T ) T 1 : (∇E) T 1 T = ∇ • ∇E T • E T -(∇E) T 3 : ∇E T . (10)
In this result, one must keep in mind that ∇E (or ∇E T ) is a third-order tensor. Since the transposition operation is not unique for this type of tensors, we have used the superscript T 1 to denote the transpose that permutes the first and second indices as defined in Eq. (B8) of Appendix B, i.e. (C is a third-order tensor), C T 1 ijk = C jik , while the superscript T 3 denotes the transpose that permutes the first and third indices as defined in Eq. (B9) of Appendix B, i.e., C T 3 ijk = C kji .

In addition, throughout this work, the nested convention is adopted for double inner products (see Eqs. (B5) and (B7) in Appendix B).

Applying the spatial averaging theorem to the first term on the rhs of the last above equation and using similar arguments as those used to simplify Eq. (8) (i.e., the nonslip-like boundary condition in Eq. (5c) and the fact that average quantities are considered as constants at the closure level), it is not hard to conclude that this term is zero. Consequently, Eq. ( 9) can now be written as

H T = (∇E) T 3 : ∇E + ρ µ E T • ∇ • v β E . ( 11
)
The first term on the rhs of this result is clearly symmetric. As a matter of fact, under creeping-flow conditions, this is the only term remaining in the definition of the intrinsic permeability tensor and it is the proof that this tensor is symmetric [START_REF] Auriault | Dynamics of porous saturated media, checking of the generalized law of Darcy[END_REF] (see also Problems 4 and 5 in the work of Whitaker 1 ). Therefore, the rest of this section is dedicated to the analysis of the second term on the rhs of Eq. ( 11) (which is originated from inertial effects at the pore-scale).

Let us commence the analysis by noting that, due to the continuity equation, the last term in Eq. ( 11) may also be written as (see Eq. ( B19))

E T • ∇ • v β E = E T v β : ∇E . ( 12
)
At this point, it is convenient to note that

E T • ∇ • v β E T = ∇ • v β E T • E = ∇ • v β E T • E -E T v β : ∇E . ( 13
)
In order to obtain the last equality, we took into account Eq. (B20). As shown above, the use of the spatial averaging theorem together with the interfacial boundary condition leads to conclude that the first term on the rhs of the above equation is null. In this way, substitution of Eq. ( 12) on the remaining term on the rhs of Eq. ( 13) yields

E T • ∇ • v β E = -E T • ∇ • v β E T . ( 14
)
This is the proof that the second term in the decomposition of the apparent permeability tensor, given in Eq. ( 11), is skewsymmetric and the tensor H can finally be written as

H = (∇E) T 3 : ∇E Symmetric part - ρ µ E T • ∇ • v β E Skew-symmetric part . ( 15
)
In this way, it can be deduced that the first term on the rhs of this last expression, which results from dissipative-like transport, is the symmetric part of the apparent permeability tensor and that the second term, arising from convective transport, is its skewsymmetric part. In other words, the decomposition shown in Eq. ( 15) readily provides the irreducible parts of the tensor H.

B. Slip flow

From the results presented in the previous paragraphs, it may be inferred that the permeability tensor remains symmetric for creeping flow, whatever the conditions. However, this conclusion should be carefully considered because it is only applicable whenever it is safe to impose nonslip conditions at the solid-fluid interface. To prove this point, consider the steady, non-inertial, Newtonian, and slightly compressible flow of a barotropic fluid saturating the pores of a rigid and homogeneous porous medium. The situation is again the one schematically represented in Fig. 1. Under isothermal conditions, the governing total mass, momentum, and state equations at the pore-scale are

∂ ρ ∂t + ∇ • (ρv β ) = 0, in V β , (16a)
0 = -∇p β + µ∇ 2 v β , in V β , (16b)
ρ = f p β , in V β . ( 16c
)
At the solid-fluid interface, the following slip-condition is applicable: 12,32-35

v β = - 2 -σ σ ξ λ β I -nn • n • ∇v β + ∇v β T , (16d)
with σ and λ β being the accommodation coefficient and the mean free path, respectively. Upscaling of this problem with the use of the volume averaging method was addressed earlier by Lasseux et al. [START_REF] Lasseux | An improved macroscale model for gas slip flow in porous media[END_REF] where the following set of macroscale equations was obtained:

∂ ρ β ∂t + ∇ • ρ β v β β = 0, (17a)
v β = - K s µ • ∇ p β β , (17b)
ρ β = f p β β , (17c)
where K s is an apparent slip-corrected permeability tensor, which can be computed from

K s = D . (18)
Here, the closure variable D is obtained from the solution of the following boundary-value problem:

∇ • D = 0, in V β , ( 19a
) 0 = -∇d + ∇ 2 D + I, in V β , (19b)
D = -ξ λ β (I -nn)• n• ∇D + (∇D) T 1 , at A βσ , (19c) ψ(r) = ψ(r + l i ); ψ = D, d, (19d)
d β = 0, (19e)
where, in the boundary condition in Eq. (19c), the superscript T 1 denotes the transpose of a third-order tensor as defined in Eq. (B8) in Appendix B. The definitions of the closure variables D and d are analogous to those of E and e presented above, respectively. In addition, λ β represents the mean free path at the intrinsic-average density and it can thus be regarded as a constant at the closure level. [START_REF] Lasseux | A macroscopic model for slightly compressible gas slip-flow in homogeneous porous media[END_REF] In Appendix B of the work by Lasseux et al., [START_REF] Lasseux | An improved macroscale model for gas slip flow in porous media[END_REF] an analysis similar to the one used in Sec. II A for the properties of the apparent permeability tensor for inertial momentum transport is provided. Here, only the final reformulation of the apparent slip-corrected permeability reported in this reference is recalled. It is given by

K s = (∇D) T 3 : ∇D - 1 V A β σ (∇D) T 3 : nD dA. ( 20
)
Just like for flow with inertia, the first term on the rhs of this expression results from the dissipation-like term in Eq. (19b) and it is not difficult to demonstrate that it is symmetric. The remaining term is, in general, non-symmetric. This observation is consistent with the study by Skjetne and Auriault 8 using the homogenization method. As a matter of fact, the last term in Eq. ( 20) is only symmetric if the solid-fluid interface exhibits symmetries across the three planes parallel to the edges of the periodic unit cell and passing through its centroid. Furthermore, after performing an order of magnitude analysis, Lasseux et al. [START_REF] Lasseux | An improved macroscale model for gas slip flow in porous media[END_REF] concluded that, for conditions in which ξKn ε (with Kn = λ β / β being the macroscale Knudsen number), K s can be considered to be quasi-symmetric. This analysis thus shows that, under creeping slip-flow conditions, it cannot be concluded that the apparent permeability tensor in Darcy's law is a symmetric tensor.

Finally, using the irreducible decomposition of any second-order tensor and the identity given in Eq. (B17), Eq. ( 20) may be expressed as

K s = (∇D) T 3 : ∇D - 1 2V A β σ (∇D) T 3 : nD + ∇D : (nD) T 3 dA Symmetric part - 1 2V A β σ (∇D) T 3 : nD -∇D : (nD) T 3 dA Skew-symmetric part , ( 21
)
where the symmetric and skew-symmetric terms are clearly identified. The existence of the skew-symmetric term may contribute to a misalignment of the velocity with the macroscopic pressure gradient as it is the case, in general, for inertial flow.

C. Slightly compressible slip-flow in a fracture

The conclusion reached in Sec. II B regarding slip flow in the bulk of a porous medium raises the question on whether it remains valid in the case of slip-flow within the confined space between two surfaces, i.e., in a fracture. The purpose of the present section is hence to investigate the properties of the effective-medium coefficient present in the upscaled model for the steady Newtonian flow of a fluid ( β-phase) saturating a fracture. The study is directed to slightly compressible slip-flow under isothermal conditions. The fracture, of local aperture h, consists of two impermeable rough walls (σ-phase) as sketched in Fig. 2(a). Contact spots may also occur between the upper and lower surfaces of the fracture as it would be the case, for instance, for the contact under normal stress between two manufactured solid surfaces (see Fig. 2(a)). When this contact exists, their contours shall be denoted by C βσ in the mean plane of the fracture (see Fig. 2

(b)). The incompressible flow version of this problem has been the object of many studies

for the leak-rate determination through mechanical seals [START_REF] Marie | An integrated approach to characterize liquid leakage through metal contact seal[END_REF][START_REF] Marie | Experimental leak-rate measurement through a static metal seal[END_REF][START_REF] Vallet | Real versus synthesized fractal surfaces: Contact mechanics and transport properties[END_REF][START_REF] Pérez-Ràfols | Modelling of leakage on metal-to-metal seals[END_REF] or for the lubrication of rough surfaces. [START_REF] Prat | Averaged Reynolds equation for flows between rough surfaces in sliding motion[END_REF] The boundary-value problem governing flow within the fracture is the one provided by Equation ( 16), except a steady version of the mass conservation Equation (16a) is considered. When the local slope of the fracture walls' roughness (i.e., tan α, see Fig. 2(a)) remains small compared to unity, the form of the above-mentioned mass and momentum conservation equations, pre-integrated in the direction of the thickness of the fracture, can be easily obtained, reducing the dimension of the problem from three to two:

∇ • q β = 0, in A βM , (22a)
q β = -ρ k µ ∇p β , in A βM , (22b)
ρ = f p β , in A βM , (22c)
q β • n = 0, at C βσ , ( 22d
)
where q β represents the mass flow rate per unit length of the fracture in each direction, while k stands for the so-called local transmissivity of the fracture aperture field that includes slip effects, i.e., Moreover, in this last equation, as in Sec. II B, λ β denotes the mean-free path of the fluid molecules at the pressure and temperature under consideration. In addition, in Eqs. (22), A βM represents the surface where no effective contact takes place and n is the unit normal vector at C βσ directed from the β-phase to the σ-phase.

k = 1 12 h 3 1 + 6 ξ λ β h . ( 23
Equations (22), which basically correspond to the Reynolds model including slip, are still operating at the scale of roughnesses, and, from a practical point of view, it is useful to provide a description of the flow over a representative elementary surface A in which the portion that excludes the effective contact area is denoted by A β . To this goal, an upscaling can be carried out using the following definitions of the averages, analog to the ones employed for volume averaging, namely, the superficial area averaging operator

ψ β s = 1 A A β ψ β dA (24a)
and intrinsic averaging operator

ψ β β s = 1 A β A β ψ β dA, ( 24b
)
where A and A β represent the areas of A and A β , respectively. Similarly, the surface version of the spatial averaging theorem takes the form

∇ψ β s = ∇ ψ β s + 1 A C β σ nψ β d . ( 25
)
When the superficial average of the problem in Eqs. ( 22) is carried out, it can be shown that the steady macroscopic Reynolds model for slightly compressible slip flow may be written as

∇ • q β s = 0, ( 26a
)
q β s = -ρ β s ε s K µ • ∇ p β β s , (26b)
ρ β s = f p β β s (26c) in which ε s = A β A = 1 -δ s ,
where δ s stands for the so-called load-bearing capacity. At this point, it is opportune to clarify that the slightly compressible hypothesis is supported by the constraint

ρ = ρ -ρ β s ρ β s . ( 27
)
In the expression of the macroscale flow-rate in Eq. (26b), K is the effective transmissivity tensor characterizing the fracture that is given by

K = k ∇b + I β s , ( 28
)
where k is defined as

k = 1 12 h 3 1 + 6 ξ λ βs h , (29)
λ βs being now the mean-free path at the average density ρ β s . In the expression of K above, b is the closure variable that maps ∇ p β β s onto the spatial variations of the fluid pressure and is the solution of the following boundary-value problem:

∇ • k ∇b + I = 0, in A β , (30a) -n • ∇b + I = 0, at C βσ , (30b) b β s = 0, (30c) b(r + l i ) = b(r), i = 1, 2. ( 30d
)
The interest shall now be dedicated to the analysis of the symmetry properties of the transmissivity tensor K. With this purpose in mind, the outer product of Eq. (30a) with b may be formed, and while taking the superficial area average of the result, one obtains

∇ • k ∇b + I b s + 1 A C β σ kn • ∇b + I bd -k ∇b T + I • ∇b s = 0. (31)
To arrive at this expression, the identity given in Eq. (B18) of Appendix B was employed as well as the averaging theorem in Eq. ( 25).

Because of periodicity and due to the slightly compressible flow assumption, the average k ∇b + I b s can be treated as a constant at the closure level so that the first term on the left-hand side of Eq. ( 31) is zero. The second term vanishes also when the boundary condition given in Eq. ( 30b) is taken into account and this yields

k∇b β s = -k∇b T • ∇b β s (32)
where we have employed the intrinsic average instead of the superficial average. When this last result is introduced back into Eq. ( 28), the effective transmissivity tensor takes the form (33) and this final result proves that K is a symmetric tensor.

K = kI β s -k∇b T • ∇b β s
As important remarks following the development performed in this section, it should be noted that:

• The expression for the effective transmissivity tensor, given in Eq. ( 33), remains completely valid if the flow takes place without slip effects and/or if the flow is perfectly incompressible. • More generally, the conclusions may be extended to the effective diffusivity tensor resulting from the upscaling of a diffusion process within a heterogeneous medium when the local diffusivity can be considered as a continuous function of space. In particular, it is readily applicable to the permeability tensor resulting from the large-scale averaging process carried out on the onephase Darcy flow in a heterogeneous porous medium when the permeability field is considered as a continuous function of space as reported by Quintard and Whitaker. 41

D. One-phase flow in a medium embedding a porous matrix and a clear fluid region

Up to this point of the investigation, momentum transport has not been considered in a heterogeneous porous medium featuring discontinuous properties and the question now remains on whether the conclusions reached so far are still applicable in this situation. To illustrate this, we are hence interested to study one-phase creeping flow of a fluid-phase β within a heterogeneous system in which a porous matrix (ωregion) and a clear fluid (η-region) can be distinctly identified as for example in a vugular (see Fig. 3) or a fractured porous medium. This problem has been analyzed by Arbogast and Lehr [START_REF] Arbogast | Homogenization of a Darcy-Stokes system modeling vuggy pous media[END_REF] and more recently by Golfier et al. [START_REF] Golfier | Investigation of the effective permeability of vuggy or fractured porous media from a Darcy-Brinkman approach[END_REF] The macroscale model used in the latter work to describe the incompressible Newtonian flow in such a material is the result of a first upscaling of the continuity and Stokes equations at the microscale and it is given by

∇ • V ω β = 0 in the ω-region, (34a)
0 = -∇P ω β + ρg + µ ε ∇ 2 V ω β -µK -1 ω • V ω β in the ω-region, (34b)
V ω β = V η β at A ωη , (34c)
n • P η β -P ω β I + µ ε -1 ∇V ω β -∇V η β + ε -1 ∇V ω β -∇V η β T = 0 at A ωη , (34d) ∇ • V η β = 0 in the η-region, (34e)
0 = -∇P η β + ρg + µ∇ 2 V η β in the η-region. (34f)
In these equations, V κ β and P κ β (κ = ω, η) are the macroscale velocity and pressure in the κ-region, ε and K ω are the porosity and intrinsic permeability tensor of the ω-porous region. In the boundary condition given in Eq. (34d), n denotes the unit normal vector at the ωη boundary A ωη , directed from the ωregion toward the η-region (i.e., n ωη as shown in Fig. 3). After applying the averaging procedure to the above equations, [START_REF] Golfier | Investigation of the effective permeability of vuggy or fractured porous media from a Darcy-Brinkman approach[END_REF] the megascale model operating on the effective medium reads

∇ • V * β = 0, (35a)
FIG. 3. Sketch of a megascopic vugular medium saturated with a single fluid phase that contains a homogeneous porous medium (ω-region) and clear fluid zones (η-region).

V * β = -

K * µ • ∇P * β -ρg , (35b)
where V * β and P * β represent the megascale average velocity and pressure defined as

Ψ * β = {Ψ β } = 1 V V Ψ β dV = {Ψ ω β } ω + {Ψ η β } η = 1 V V ω Ψ ω β dV + 1 V V η Ψ η β dV , (36)
where V κ (κ = ω, η) represents the portion of the κ-region contained within the averaging volume V = V ω + V η . In the megascale momentum Equation (35b), K * represents the effective permeability tensor that can be computed from the solution of the associated closure problem given by

∇ • D ω = 0, in V ω , ( 37a
)
0 = -∇d ω + 1 ε ∇ 2 D ω -K -1 ω • D ω + I, in V ω , (37b) B.C. 1 D ω = D η at A ωη , (37c)
B.C. 2 n • -Id ω + 1 ε ∇D ω + ∇D ω T 1 = n • -Id η + ∇D η + ∇D η T 1 at A ωη , (37d)
∇ • D η = 0, in V η , ( 37e
)
0 = -∇d η + ∇ 2 D η + I, in V η , ( 37f
)
d * = 0, (37g)
D * = K * , ( 37h
)
ψ k (r + l i ) = ψ k (r), ψ = d, D, k = ω, η, i = 1, 2, 3. (37i)
In these equations, D κ and d κ (κ = ω, η) are closure variables that map V * β onto the velocity and pressure variations in the κ-region, respectively, whereas D * and d * are the megascale averages based on the definition given in Eq. (36).

The analysis of the symmetry properties of the effective permeability tensor K * starts by first reformulating the momentum-like Equations (37b) and (37f), which, while taking into account the fact that D κ , (κ = ω, η) are solenoidal tensor fields as indicated by Eqs. (37a) and (37e), may indeed be written as

0 = -∇d ω + 1 ε ∇ • ∇D ω + ∇D ω T 1 -K -1 ω • D ω + I (38)
and

0 = -∇d η + 1 ε ∇ • ∇D η + ∇D η T 1 + I. ( 39
)
When the former of these two equations is pre-multiplied by D T ω and when the ω-regional average of the result is formed, the following expression is obtained:

0 = 1 V A ω η D T ω • nd ω dA + 1 ε 1 V × A ω η n • ∇D T ω + ∇D T ω T 3 • D ω T dA - 1 ε ∇D T ω + ∇D T ω T 3 T 1 : ∇D ω T 1 T ω -D T ω • K -1 ω • D ω ω + D T ω ω . ( 40
)
This last result is subject to the use of the averaging theorem for the ω-regional average, which has exactly the same form as the one given in Eq. (A2), together with the fact that d ω and D ω are periodic. Moreover, the first term on the rhs of the above equation results from the fact that D ω is solenoidal. The second term on the rhs of Eq. ( 40) results from the use of the identities given in Eqs. (B12) and (B14) in which the arbitrary twice-derivable second-order tensor A is identified as D ω . It also makes use of the identity given in Eq. (B22), in which ∇A is identified as

∇D T ω + ∇D T ω T 3
and B to D ω .

At this point, it is worth noticing that n

• ∇D T ω T = n • ∇D ω and n • ∇D T ω T 3 T = n • ∇D ω T 1
, so that taking into account Eqs. (B15) and (B16), Eq. ( 40) can be rewritten as

0 = 1 V A ω η D T ω • n • -Id ω + 1 ε ∇D ω + ∇D ω T 1 dA - 1 ε ∇D ω T 3 : ∇D ω + ∇D ω T 1 T ω -D T ω • K -1 ω • D ω ω + D T ω ω . ( 41
)
Let now attention be directed toward the η-region. Using the same procedure employed above, it follows that

0 = - 1 V A η η D T η • n • -Id η + ∇D η + ∇D η T 1 dA -∇D η T 3 : ∇D η + ∇D η T 1 T η + D T η η . ( 42
)
When Eqs. ( 41) and ( 42

K * = 1 ε ∇D ω T 3 : ∇D ω + ∇D ω T 1 ω + ∇D η T 3 : ∇D η + ∇D η T 1 η + D T ω • K -1 ω • D ω ω . ( 43
)
In this last expression, we have made use of the fact that K ω is a symmetric tensor to infer that D T ω • K -1 ω • D ω ω is also symmetric. To reach this conclusion, the symmetry property of K ω , which is known in the creeping flow regime with the no-slip boundary condition at the fluid-solid interface, as mentioned earlier, is indeed fundamental. This is a clear example showing why the knowledge of the symmetry properties of transport coefficients at the macroscale is extremely important. Clearly, the first two terms in Eq. ( 43) are also symmetric as shown in Appendix B (see Eqs. (B25) and (B26)) making the effective permeability tensor K * a symmetric tensor and this represents the final important result of this section.

E. Creeping two-phase flow in homogeneous porous media

So far, the analysis of macroscopic coefficients has been restricted to single phase-flow in porous structures. As a final momentum transport case, let us then consider incompressible, immiscible, and Newtonian two-phase momentum transport in the bulk region of a rigid porous medium as sketched in Fig. 4. The governing pore-scale equations and interfacial boundary conditions for this transport process are given by

∇ • v k = 0, in the k-phase (k = β, γ), (44a)
0 = -∇p k + µ k ∇ 2 v k , in the k-phase(k = β, γ), (44b) B.C. 1 v k = 0, at A kσ (k = β, γ), (44c)
B.C. 2 v β = v γ , at A βγ , (44d)
B.C. 3 n βγ • T β = n βγ • T γ + 2σHn βγ , at A βγ . (44e)
Here, the nonslip condition has been imposed for both fluid phases in contact with the solid phase. In addition, assuming that no mass transport is taking place between the fluid phases, it is reasonable to impose continuity of the velocity fields at A βγ . Moreover, at this interface, the normal stress is exactly compensated by the Laplace pressure as expressed in Eq. (44e) where the total stress tensor is denoted as T k = -Ip k + µ k ∇v k + (∇v k ) T (k = β, γ), while σ and H represent the interfacial tension and the mean curvature of A βγ , respectively. This problem has been studied by Whitaker [START_REF] Whitaker | The closure problem for two-phase flow in homogeneous porous media[END_REF][START_REF] Whitaker | Flow in porous media II: The governing equations for immiscible, two-phase flow[END_REF] and by Lasseux et al. [START_REF] Lasseux | Determination of permeability tensors for two-phase flow in homogeneous porous media: Theory[END_REF] using the volume averaging method. The resulting upscaled model for momentum transport can be expressed as follows: This model was also derived using the homogenization method by Auriault. [START_REF] Auriault | Nonsaturated deformable porous media: Quasistatics[END_REF] In Eqs. (45), the permeability tensors K * αk (α, k = β, σ) are defined in terms of the corresponding closure variables as

v β = - K * ββ µ β • ∇ p β β - K * βγ µ γ • ∇ p γ γ , ( 45a
)
v γ = - K * γγ µ γ • ∇ p γ γ - K * γβ µ β • ∇ p β β . ( 45b
K * αk = -D αk , α, k = β, γ. (46)
Following the work of Lasseux et al., [START_REF] Lasseux | Determination of permeability tensors for two-phase flow in homogeneous porous media: Theory[END_REF] we shall refer to the tensors K * ββ and K * γγ as the dominant permeability tensors and to the tensors K * βγ and K * γβ as the coupling permeability tensors. The closure variables D ββ and D γβ are the solution of the following boundary-value problem.

Problem I

∇ • D kβ = 0, in V k (k = β, γ), (47a)
-∇d kβ + ∇ • ∇D kβ + ∇D kβ T 1 = δ kβ I, in V k (k = β, γ), (47b)
B.C. 1 D kβ = 0, at A kσ (k = β, γ), (47c)
B.C. 2 D ββ = D γβ , at A βγ , (47d)
B.C. 3 µ β n βγ • -Id ββ + ∇D ββ + ∇D ββ T 1 = µ γ n βγ • -Id γβ + ∇D γβ + ∇D γβ T 1 , at A βγ , (47e)
ψ kβ (r + l i ) = ψ kβ (r), i = 1, 2, 3; ψ = d, D; k = β, γ, (47f)
d kβ k = 0, k = β, γ. (47g)
In Eq. (47b), δ kβ is the Dirac delta function. Note that in this equation, we have taken into account the solenoidal nature of the tensors

D kβ (k = β, γ), so that ∇ ∇ • D kβ = ∇ • ∇D kβ T 1 = 0.
In a similar manner, the closure variables D βγ and D γγ must satisfy the following closure problem.

Problem II

∇ • D kγ = 0, in V k (k = β, γ), (48a)
-∇d kγ + ∇ • ∇D kγ + ∇D kγ T 1 = δ kγ I, in V k (k = β, γ), (48b)
B.C. 1 D kγ = 0, at A kσ (k = β, γ), (48c)
B.C. 2 D βγ = D γγ , at A βγ , (48d)
B.C. 3 µ β n βγ • -Id βγ + ∇D βγ + ∇D βγ T 1 = µ γ n βγ • -Id γγ + ∇D γγ + ∇D γγ T 1
, at A βγ , (48e)

ψ kγ (r + l i ) = ψ kγ (r), i = 1, 2, 3; ψ = d, D; k = β, γ, (48f)
d kγ k = 0, k = β, γ. (48g)
These closure problems are restricted to cases in which the surface tension has no effects at the closure problem level. [START_REF] Auriault | Nonsaturated deformable porous media: Quasistatics[END_REF][START_REF] Whitaker | Flow in porous media II: The governing equations for immiscible, two-phase flow[END_REF] This is a reasonable assumption whenever the capillary number is much smaller than unity; i.e., Ca = µ /σ 1, with µ being the largest of O(µ β v β β ) or O(µ γ v γ γ ). To commence the derivations, let attention be focused to Problem I, in specific, consider Eq. (47b) for k = β, premultiply it by D T ββ , and apply the superficial average defined in Eq. (A1a) of Appendix A on the resulting equation in order to obtain

-D T ββ • ∇d ββ β + D T ββ • ∇ • ∇D ββ + ∇D ββ T 1 β + K * T ββ = 0. (49)
Here, we have used the subscript β in the superficial average to clearly indicate that this average is taken only in the β-phase.

At this point, one may follow the same line of derivations performed in Section II D in order to obtain

0 = 1 V A βγ D T ββ • n βγ • -Id ββ + ∇D ββ + ∇D ββ T 1 dA -∇D ββ T 3 : ∇D ββ + ∇D ββ T 1 T β + K * T ββ . (50a)
Let us return the attention to Eq. (47b) but now for k = γ with a pre-multiplication by D T γβ ; performing the same procedure used just above, it is not hard to deduce that

0 = - 1 V A βγ D T γβ • n βγ • -Id γβ + ∇D γβ + ∇D γβ T 1 dA -∇D γβ T 3 : ∇D γβ + ∇D γβ T 1 T γ . ( 50b
)
When this last result is multiplied by µ γ and it is added to Eq. (50a) multiplied by µ β , the resulting equation takes the form

K * ββ = µ γ µ β ∇D γβ T 3 : ∇D γβ + ∇D γβ T 1 γ + ∇D ββ T 3 : ∇D ββ + ∇D ββ T 1 β . (51)
Here, the interfacial boundary conditions given in Eqs. (47d) and (47e) were taken into account. The two average terms on the rhs of this equation can be easily shown to be symmetric (see the proof in Appendix B, Eqs. (B25) and (B26)), thus making K * ββ a symmetric tensor. This conclusion is consistent with those obtained separately by Auriault 9,15 and Whitaker. [START_REF] Whitaker | The closure problem for two-phase flow in homogeneous porous media[END_REF] The above procedure can be repeated with Eq. (47b) (for k = β) when it is pre-multiplied by D T βγ and also with this same equation, when evaluated for k = γ and pre-multiplied by D T γγ to obtain, after few steps, which are analogues to those used above,

K * βγ = µ γ µ β ∇D γβ + ∇D γβ T 1 T 3 : ∇D γγ T 1 γ + ∇D ββ + ∇D ββ T 1 T 3 : ∇D βγ T 1 β . ( 52
)
To arrive at this result, we have made use of the identity given in Eq. (B13).

The remaining tensors in the macroscale model can be obtained by performing a similar analysis to Problem II. For the sake of brevity in presentation, we only report the final results, which are given by

K * γβ = µ β µ γ ∇D βγ + ∇D βγ T 1 T 3 : ∇D ββ T 1 β + ∇D γγ + ∇D γγ T 1 T 3 : ∇D γβ T 1 γ , (53a)
K * γγ = µ β µ γ ∇D βγ T 3 : ∇D βγ + ∇D βγ T 1 β + ∇D γγ T 3 : ∇D γγ + ∇D γγ T 1 γ . ( 53b
)
This last result on K * γγ clearly exhibits the same symmetry properties as K * ββ (see Eq. ( 51)). Unfortunately, the tensors K * γβ and K * βγ do not show any particular symmetry properties. Nevertheless, a reciprocity relationship can still be derived. To prove this point, it is useful to consider the identities given in Eqs. (B23) and (B24), so that

µ β K * βγ = µ γ K * T γβ . ( 54
)
The reciprocity relationship given in Eq. ( 54), which is in agreement with Onsager's arguments, has been previously obtained by Auriault 9,[START_REF] Auriault | Nonsaturated deformable porous media: Quasistatics[END_REF] and Lasseux et al. [START_REF] Lasseux | Determination of permeability tensors for two-phase flow in homogeneous porous media: Theory[END_REF] As a summary for this case study, it is worth recalling that, for the upscaled model given in Eqs. (45), the tensors K * ββ and K * γγ are symmetric. In addition, the coupling permeability tensors K * βγ and K * γβ , although they do not exhibit a special symmetry or skew-symmetry property, are related as shown in Eq. ( 54), thus reducing the number of effective medium coefficients to be computed to 3 instead of 4. The irreducible decomposition of second-order tensors can be used for the coupling tensors in order to express them in terms of their symmetric and skew-symmetric parts. However, for the sake of brevity, this decomposition is not presented here. As a final note, it is worth pointing out that the relevance of the coupling permeability tensors has been addressed in several works, [START_REF] Zarcone | Détermination expérimentale du couplage visqueux dans les écoulements diphasiques en milieu poreux[END_REF][START_REF] Ramakrishnan | Measurement of off-diagonal transport coefficients in two-phase flow in porous media[END_REF][START_REF] Li | Pore-scale investigation of viscous coupling effects for two-phase flow in porous media[END_REF] from which a definite conclusion is still lacking.

III. DISPERSIVE HEAT AND MASS TRANSPORT

This last part of the article is devoted to the analysis of the effective transport coefficients involved in the upscaled equations for dispersive heat and mass transport in the bulk region of a rigid and homogeneous porous medium. The configuration is the one depicted in Fig. 1 showing the different scales and the domain over which the pore-scale equations are averaged. Heat transport, involving conduction and convection, is first considered under non-equilibrium conditions and subsequently passive mass dispersion is addressed.

A. Heat transport

In their study of heat dispersion in homogeneous porous media, Quintard et al. [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: Numerical results for effective properties[END_REF] considered the following pore-scale equations:

ρc p β ∂T β ∂t + ρc p β v β • ∇T β = ∇ • k β ∇T β , in V β , ( 55a
)
ρc p σ ∂T σ ∂t = ∇ • (k σ ∇T σ ) , in V σ , (55b)
where convection only takes place in the fluid-phase experiencing an incompressible and Newtonian flow governed by the Stokes equations. The above equations are coupled at the solid-fluid interface by means of the following boundary conditions:

T β = T σ , at A βσ , ( 55c
) n • k β ∇T β = n • k σ ∇T σ , at A βσ . ( 55d
)
The result of applying the volume averaging method, without the imposition of the local thermal equilibrium assumption, is the following set of upscaled equations:

ε β ρc p β ∂ T β β ∂t + ε β ρc p β v β β • ∇ T β β -u ββ • ∇ T β β -u βσ • ∇ T σ σ = ∇ • K ββ • ∇ T β β + K βσ • ∇ T σ σ -a h T β β -T σ σ , (56a)
ε σ ρc p σ ∂ T σ σ ∂t -u σβ • ∇ T β β -u σσ • ∇ T σ σ = ∇ • K σβ • ∇ T β β + K σσ • ∇ T σ σ + a h T β β -T σ σ , ( 56b
)
where the scalar a h, the vectors u αk , and the tensors K αk (α, k = β, σ) are defined in terms of the corresponding closure variables (see Eq. ( 19) in the work of Quintard et al. [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: Numerical results for effective properties[END_REF] ). In this section, the attention is focused on the effective thermal conductivity tensors defined as

K ββ = k β ε I + 1 V A β σ nb ββ dA -ρc p β ṽβ b ββ , (57a)
K βσ = k β V A β σ nb βσ dA -ρc p β ṽβ b βσ , ( 57b
)
K σβ = - k σ V A β σ nb σβ dA, ( 57c
)
K σσ = k σ         (1 -ε)I - 1 V A β σ nb σσ dA         . ( 57d
)
The closure variables b ββ and b σβ are the solution of the following boundary-value problem.

In this way, Eq. (57a) can now be written as

K ββ = k β εI -k β ∇b T ββ • ∇b ββ β -k σ ∇b T σβ • ∇b σβ σ -(ρc p ) β v β • ∇b ββ b ββ β -2(ρc p ) β ṽβ b ββ β . ( 70
)
Taking into account the spatial decomposition of the velocity and the average constraint in Eq. (58f), we may group the last two terms of the above equation into a single one and express K ββ in its final form,

K ββ = k β εI -k β ∇b T ββ • ∇b ββ β -k σ ∇b T σβ • ∇b σβ σ -(ρc p ) β v β • (∇b ββ + 2I)b ββ β . (71)
Furthermore, on the basis of the interfacial boundary condition in Eq. (58c), it can also be deduced that

k β k σ K σβ = k σ ∇b T σβ • ∇b σβ σ + k β ∇b T ββ • ∇b ββ β + (ρc p ) β v β • (∇b ββ + I)b ββ β . ( 72
)
With the aim of analyzing the remaining effectivemedium coefficients, let us now make the outer product of Eq. (60a) with b βσ and apply the superficial averaging operator for the β-phase. After performing similar steps to those used above, one easily arrives at

∇ • (k β ∇b βσ )b βσ β = k σ V A β σ nb σσ dA + 1 V A β σ n • (k σ ∇b σσ ) b σσ dA -k β ∇b T βσ • ∇b βσ β . (73)
To make further progress, let us make the outer product of Eq. (60b) with b σσ and apply the superficial average in the σ-phase. After few steps that also involve the use of Eq. ( 73), the following expression is obtained:

k σ V A β σ nb σσ dA = (ρc p ) β v β • ∇b βσ b βσ β + k β ∇b T βσ • ∇b βσ β + k σ ∇b T σσ • ∇b σσ σ . (74)
Substitution of this result into Eq. (57d) yields

K σσ = k σ (1 -ε)I -k β ∇b T βσ • ∇b βσ β -k σ ∇b T σσ • ∇b σσ σ -(ρc p ) β v β • ∇b βσ b βσ β . (75)
Finally, taking into account the interfacial boundary condition given in Eq. (60c), the result given in Eq. ( 74) can also be used to express K βσ as follows:

k σ k β K βσ = k β ∇b T βσ • ∇b βσ β + k σ ∇b T σσ • ∇b σσ σ + (ρc p ) β v β • ∇b βσ - k σ k β I b βσ β . ( 76
)
Observation of the derived expressions for all the effective thermal conductivity coefficients shows that they all contain symmetric terms and one volume-averaged term that vanishes under non-convective conditions. It is thus interesting to carry out a deeper analysis in order to determine the structure of this last term. As a matter of fact, for all the effective conductivities, this term involves the following average, v β • ∇b βk b βk β (k = β, σ), which, due to the solenoidal nature of the velocity field, can be written in the following alternative form:

v β • ∇b βk b βk β = ∇ • v β b βk b βk β = ∇ • v β b βk b βk β -b βk v β • ∇b βk β . (77)
Here we have also taken into account the identity given in Eq. (B18). It is not hard to prove that the first term on the rhs of the last equality of the above expression is null as a result of the use of the spatial averaging theorem together with the nonslip boundary condition at the β-σ interface. Finally, taking once again into account the solenoidal nature of v β , it follows that

∇ • (v β b βk )b βk β = -b βk ∇ • (v β b βk ) β = -∇ • (v β b βk )b βk T β , (78)
hence showing that this term is skew-symmetric. Unfortunately, since the outer product between the velocity and the closure variable b βk has no particular symmetry properties, the term v β b βk β , which is present in all the effective thermal conductivities (except in K σσ ) is neither symmetric nor skewsymmetric. However, this term can be simply decomposed into its irreducible form. In this way, all the thermal conductivity tensors can be expressed as

K ββ = k β εI -k β ∇b T ββ • ∇b ββ β -k σ ∇b T σβ • ∇b σβ σ -(ρc p ) β v β b ββ + b ββ v β β Symmetric part -(ρc p ) β v β • ∇b ββ b ββ β -(ρc p ) β v β b ββ -b ββ v β β Skew-symmetric part , ( 79a
)
k β k σ K σβ = k σ ∇b T σβ • ∇b σβ σ + k β ∇b T ββ • ∇b ββ β + 1 2 (ρc p ) β v β b ββ + b ββ v β β Symmetric part + ρc p β v β • ∇b ββ b ββ β + 1 2 (ρc p ) β v β b ββ -b ββ v β β Skew-symmetric part , (79b)
k σ k β K βσ = k β ∇b T βσ • ∇b βσ β + k σ ∇b T σσ • ∇b σσ ) σ - k σ 2k β (ρc p ) β v β b βσ + b βσ v β β Symmetric part + ρc p β v β • ∇b βσ b βσ β - k σ 2k β (ρc p ) β v β b βσ -b βσ v β β Skew-symmetric part , (79c)
K σσ = k σ (1 -ε)I -k β ∇b T βσ • ∇b βσ β -k σ ∇b T σσ • ∇b σσ σ Symmetric part -ρc p β v β • ∇b βσ b βσ β Skew-symmetric part . ( 79d
)
From the above results, the following comments are in order.

(1) The analysis carried out here allows us to exhibit the dependence on convection of all the coefficients. This contrasts with the original form given in Eqs. (57), where the effective coefficients associated with the σphase (K σσ and K σβ ) do not contain a hydrodynamic dispersion contribution as it is the case for the K ββ and K βσ coefficients. (2) The reformulation of the thermal conductivity tensors shows that none of them are symmetric. However, the present developments operate decomposition into the symmetric and skew-symmetric parts in each tensor. The skew-symmetric parts of the tensors only contain convective transport terms in the β-phase. (3) For transport conditions in which the thermal Péclet number, defined as Pe = (ρc P ) β v β β β /k β , is much smaller than unity, all the coefficients are quasisymmetric and they are perfectly symmetric under purely conductive conditions. (4) For situations in which the assumption of local thermal equilibrium is reasonable, there is only one effective thermal conduction coefficient, K eff = K ββ + K βσ + K σβ + K σσ (see Eq. ( 15) in the work of Quintard et al. [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: Numerical results for effective properties[END_REF] ) to which all the above observations are applicable. (5) As explained by Whitaker 1 (see Problems 3-7), only the symmetric parts of the effective medium coefficients influence the transport equation. However, if the macroscopic boundary conditions are of Neumann or Robin-type, then the skew-symmetric parts of the tensors cannot be discarded.

B. Mass dispersion

Finally, from the previous derivations, it is of interest to address the passive mass dispersion problem in the bulk region of a rigid and homogeneous porous medium. The governing equation involving diffusive and convective mass transport of species A at the pore-scale is

∂c Aβ ∂t + v β • ∇c Aβ = ∇ • D β ∇c Aβ , in V β . (80a)
Assuming the solid phase to be impermeable to mass transport, the following interfacial boundary condition is applicable:

-n • D β ∇c Aβ = 0, at A βσ . ( 80b
)
As in the heat transport analysis, this problem is coupled to momentum transport and the fluid velocity is subject to nonslip conditions at the fluid-solid interface. As detailed by Whitaker 1 (see Chapter 3), the result from using the volume averaging method to this problem is the following upscaled model:

ε ∂ c Aβ β ∂t + v β • ∇ c Aβ β = εD * β : ∇∇ c Aβ β , (81)
where the total dispersion tensor is defined as

D * β D β = I + 1 V β A β σ nf β dA - 1 D Aβ ṽβ f β β . ( 82
)
In order to predict the values of this coefficient, it is thus necessary to solve the following closure problem:

ṽβ + v β • ∇f β = D β ∇ 2 f β , in V β , (83a)
n • ∇f β = -n, at A βσ , (83b)
f β (r + l i ) = f β (r), i = 1, 2, 3, (83c)
f β β = 0. (83d)
As recently shown by Valdés-Parada et al., [START_REF] Valdés-Parada | A new formulation of the dispersion tensor in homogeneous porous media[END_REF] the dispersion tensor can be re-written in an equivalent dimensional form, using exactly the same procedure employed in the previous paragraphs for heat transfer, as follows:

D * β D β = I -(∇f β ) T • ∇f β β - 1 D β v β • (∇f β )f β β - 2 D β v β f β β . ( 84
)
The last term in this equation is exactly twice the so-called hydrodynamic dispersion tensor in the original work by Whitaker. [START_REF] Whitaker | The Method of Volume Averaging[END_REF] This term has no specific symmetry properties in the general case, except when the periodic unit cell representative of the structure is symmetric and when the flow is along one of the symmetry axes. The third term on the rhs of Eq. ( 84) can be easily proven to be skew-symmetric and this can be obtained by following exactly the same procedure and arguments as those used for Eqs. ( 77) and (78). Certainly, after using the irreducible decomposition in this last term, the total dispersion tensor may be expressed as

D * β D β = I -(∇f β ) T • ∇f β β - 1 D β v β f β + f β v β β Symmetric part - 1 D β v β • (∇f β )f β β - 1 D β v β f β -f β v β β Skew-symmetric part . ( 85
)
Here, it is not hard to identify the same properties as for the effective thermal conductivity tensors, namely, a symmetric part present in the first three terms on the rhs of Eq. ( 85) and a skew-symmetric part (fourth and fifth terms). Not surprisingly, the skew-symmetric part of D * β only contains convective terms. There has been some confusion in the literature not only about the definition of the hydrodynamic dispersion part of D * β but also about the properties of the total dispersion tensor itself. For example, in the works by Bear et al., [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Bear | Effects of material symmetry on the coefficients of transport in anisotropic porous media[END_REF] the total dispersion tensor is argued to be totally symmetric on the basis of the conjugated thermodynamic force and flux relation. The present development evidences that this cannot be the case due to the presence of skew-symmetric terms, which are only negligible under poorly or non-convective conditions. This conclusion is in agreement with the results presented by Auriault et al. [START_REF] Auriault | On the asymmetry of the dispersion tensor in porous media[END_REF] using the homogenization method and confirmed by numerical simulations. As remarked by these authors, Onsager's relations cannot be invoked here due to the irreversible character of the mass and momentum transport processes taking place at the pore-scale.

It should be noted that the symmetric part of D * β in Eq. ( 85) is the one resulting from the method of moments, as reported by Brenner [START_REF] Brenner | Dispersion resulting from flow through spatially periodic porous media[END_REF] (see Sections 5 and 6 therein) and later on by Salles et al. [START_REF] Salles | Taylor dispersion in porous media. Determination of the dispersion tensor[END_REF] This is consistent with what has been pointed out by Koch and Brady, [START_REF] Koch | The symmetry properties of the effective diffusivity tensor in anisotropic porous media[END_REF] who showed that the method of moments, or the Lagrangian approach, can only predict the symmetric part of the total dispersion tensor. Importantly, the present analysis clearly evidences the complementary (i.e., the skew-symmetric) parts of this tensor that are missing from the method of moments analysis. As a matter of fact, these parts of D * β do not contribute to the macroscopic mass transport as can be easily inferred from the rhs of Eq. (81); this is simply due to the fact that ∇∇ c Aβ β is a symmetric tensor.

Nevertheless, as pointed out during the analysis of heat transport, this result is not extensible to situations in which mass transport near the porous medium boundaries is prescribed by Neumann or Robin-type boundary conditions. [START_REF] Whitaker | The Method of Volume Averaging[END_REF] Indeed, in this situation the skew-symmetric part of the dispersion tensor cannot be neglected in the estimation of mass flux and this is even the dominant part under strongly convective conditions (i.e., for sufficiently large values of the mass Péclet number, say Pe = v β β β /D β ). In this way, the conclusions from the work by Bear et al. [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF][START_REF] Bear | Effects of material symmetry on the coefficients of transport in anisotropic porous media[END_REF] can only hold if their total dispersion tensor is meant to be the same as the one arising from the method of moments.

IV. CONCLUSION

In this work, we have used the framework of volume averaging as an upscaling tool to study the symmetry properties of effective-medium coefficients arising in macroscopic equations in many transport processes. In all cases we have followed a similar procedure, which can be summarized as follows.

(1) Use the governing equations for the closure variables that define the macroscopic transport coefficients and perform an appropriate (inner or outer) product with the adequate closure variable. Since the definition of the effective-medium coefficient is to be recovered from the differential equation that governs the fields of the closure variable, the appropriate product to be used in each case is the one that leads to a tensorial equation of the same rank as the desired coefficient. For example, in momentum transport the corresponding differential equations are already tensorial; therefore, the corresponding product is an inner product, whereas for heat and mass transport the equation is vectorial and thus an outer product is in order. (2) Apply the superficial average to the previous result along with the averaging theorem and corresponding boundary and periodicity conditions in the closure problem to obtain, after few algebraic manipulations involving tensor analysis, an alternative expression of the corresponding macroscopic transport coefficient in terms of the closure variables. (3) The resulting expression for the effective-medium coefficients has the nice feature that it allows a straightforward identification of the symmetry and skewsymmetric parts of the tensors.

The technique outlined above was widely illustrated over several study cases, from which the following comments arise.

• For inertial one-phase flow in homogeneous porous media, the method performs the irreducible decomposition of the apparent permeability tensor, showing that the symmetric part results from viscous dissipation, while the skew-symmetric part originates from inertial transport only. Hence, under creeping flow conditions, the apparent permeability tensor corresponds to the intrinsic one, which is fully symmetric as it is well-known in the literature.

• When the slip-boundary condition is applicable, even in the creeping-flow regime, the apparent permeability tensor is not symmetric. • For one-phase slightly compressible gas flow, in the slip regime within a fracture, the effective transmissivity tensor is symmetric despite the presence of slip effects. The result obtained in this case can be generalized to any diffusion process involving continuously position dependent diffusivities. • The effective permeability tensor for one-phase flow in a heterogeneous porous medium embedding a porous matrix and a clear fluid region was shown to be symmetric when the flow is assumed to be governed by the Darcy-Brinkman equations in the porous medium and by the Stokes equations in the clear fluid region. • In the case of two-phase creeping flow in homogeneous porous media, the dominant permeability tensors are shown to be symmetric, while the coupling permeability tensors do not exhibit specific symmetry properties. A reciprocity relationship between the latter was recovered in consistency with a previous analysis. • While studying conductive and convective heat transfer under non-local thermal equilibrium conditions, it was demonstrated that the effective thermal conductivity tensors include a symmetric and a skew-symmetric part, the latter involving convective effects only. When convection is unimportant or absent, the effective thermal conductivity tensors become symmetric. • The results from the previous item are extensible to passive mass dispersion in homogeneous porous media. More importantly, our developments allow us to explicitly identify the skew-symmetric part of the total dispersion tensor that is otherwise missing when derived from the method of moments.

Beyond the study cases included in this work, the technique used here can be applied to any other transport process that is upscaled using the volume averaging method. Hence, the procedure can be regarded to be systematic for the analysis of the symmetry properties of macroscopic transport coefficients as it was the goal of this work. Finally, it is worth mentioning that the conclusions reached in this analysis remain valid whatever the microstructure of the porous medium. Indeed, the upscaling process does not require any specification of the porous medium geometry except that, at the closure level, the structure is made pseudo-periodic.

ACKNOWLEDGMENTS

Inspiration of this development is intimately related to wonderful exchanges with F.T.C.P. to whom D.L. would like to warmly dedicate this work. F.J.V.-P. is thankful to the CNRS for the facilities provided to carry out a scientific stay with D.L. during spring 2016. D.L. acknowledges the financial support from CONACyT to perform a two-week stay in UAM-Iztapalapa to work with F.J.V.-P. during autumn 2016.

APPENDIX A: OUTLINE OF THE VOLUME AVERAGING METHOD

In this section, the method of volume averaging is briefly outlined. The steps presented here are consistent with the monograph by Whitaker [START_REF] Whitaker | The Method of Volume Averaging[END_REF] and also with the more recent version reported by Wood and Valdés-Parada. [START_REF] Wood | Volume averaging: Local and nonlocal closures using a Green's function approach[END_REF] In Fig. 5, the main steps, assumptions, and tools involved in this upscaling method are schematized and listed below.

(1) First, a set of starting assumptions are introduced with the aim of defining the transport problem at the pore scale. (2) An averaging domain V (of measure V ) that contains portions of all the phases involved in the system is then defined along with the corresponding superficial

ψ α = 1 V V α ψ α dV (A1a)
and intrinsic

ψ α α = 1 V α V α ψ α dV (A1b)
averaging operators for a piece-wise continuous function, ψ α , defined everywhere in the α-phase. Actually, these averaging operators are coupled by means of the Dupuit-Forchheimer relation ψ α = ε α ψ α α , with ε α = V α /V being the volume fraction of the α-phase within the averaging domain. For the case in FIG. 5. Scheme of the upscaling process using the method of volume averaging.

which there is only one fluid phase (the β-phase) (see Fig. 1) saturating the porous matrix, the fluid volume fraction corresponds to the porosity, and throughout this work, it is denoted by ε β ≡ ε.

(3) The superficial averaging operator defined in Eq. ( A1a) is applied to the pore-scale equations. In addition, the interchange of spatial integration and differentiation can be achieved by means of the spatial averaging theorem, [START_REF] Howes | The spatial averaging theorem revisited[END_REF] which, for any quantity ψ α defined in the α-phase, is given by

∇ψ α = ∇ ψ α + 1 V A α n α ψ α dA (A2)
and a completely similar expression in its divergence form.

It is worth adding that the application of the averaging theorems often allows substitution of the corresponding interfacial boundary conditions. In Eq. (A2), A α denotes all the surfaces with which the α-phase is in contact and n α is a unit normal vector directed from the α-phase toward each of the other phases in contact. For the particular case in which there is a single fluid-phase saturating the porous medium, A β = A βσ and n βσ ≡ n. (4) At this point, the average equations are expressed in terms of both average and pointwise quantities. To eliminate the latter, the spatial decomposition introduced by Gray [START_REF] Gray | A derivation of the equations for multiphase transport[END_REF] is used,

ψ α = ψ α α + ψα , ( A3
)
ψα being the spatial deviations of ψ α about its intrinsic average. The resulting expression is the unclosed average model because it lacks of a relation between the spatial deviations and average quantities. Following the work of Whitaker, 1 the (surface and/or volumetric) integrals containing deviation terms may be conceived as filters of information coming from the pore scale. (5) With the aim of closing the average model, a boundaryvalue problem for the spatial deviations is derived and formally solved. Since no assumptions have been imposed so far, let us refer to this boundary-value problem as the exact closure problem. However, the price to be paid is that its complexity is the same as (or even greater than) that of the pore-scale model. It is thus convenient to introduce a first set of scaling postulates with the aim of filtering out the redundant information present in the pore-scale. In typical volume averaging applications, these scaling postulates consist of a set of reasonable length-scale constraints and assumptions among which it is worth recalling the following one:

r 0 L (A4)
in which r 0 is the characteristic size of the averaging domain, while and L represent the largest characteristic length of the pore-scale and the smallest characteristic length of the macroscale, respectively.

Using orders of magnitude estimates and taking into account this separation of length scales, the differential equations and boundary conditions in the closure problem can be considerably simplified. In this way, the average quantities present in the closure problem are position-invariant. A corollary of this approximation is that the deviations fields must satisfy the following average constraint:

ψα α = 0. (A5)
A convenient, although non-mandatory, simplification of the closure problem is that, in many situations, its solution domain can be reduced to a single periodic unit cell. Therefore, one may impose the following boundary condition at the entrances and exits of the unit cell: ψα (r + l i ) = ψα (r), i = 1, 2, 3, (A6)

with r and l i being a position vector and each of the lattice vectors in the unit cell, respectively. Under these conditions, let us refer to this version of the closure problem as the simplified closure problem, which has the nice quality of requiring knowledge of less information than the pore-scale model. ( 6) Substitution of the formal solution of the simplified closure problem into the filters terms of the unclosed average model leads to its closed version, which can be subsequently reduced by performing orders of magnitude analyses taking into account the separation of length scales already imposed in the simplified closure problem. After performing this final development, the resulting model may be referred to as the simplified upscaled model.

)

  FIG.1. Sketch of a system, consisting of a porous medium saturated with a single fluid phase including characteristic lengths and a sample of the averaging domain.

  FIG. 2. (a) Sketch of a fracture between two rough surfaces for which the aperture field is h. The local slope is such that tan α 1 everywhere. (b) Top view of the fracture with contact areas which contours are C β σ and averaging surface A.

  ) are added and when the boundary conditions B.C.1 and B.C.2 in Eqs. (37c) and (37d) are employed together with the definition of the effective permeability tensor in Eq. (37h), the following expression for K * arises:

)FIG. 4 .

 4 FIG.4. Sketch of a porous medium saturated with two immiscible fluid phases.

Problem I

In the above expressions,

In addition, the closure variables b βσ and b σσ are given by the solution of the following boundary-value problem.

Problem II

In this case,

Note that, at the closure problem level, all the c αk (α, k = β, σ) vectors are constants within the unit cell. Let us commence by analyzing the properties of the tensor K ββ . From Eq. (57a), it follows that, in order to study the nature of this tensor, it is convenient to make the outer product of Eq. (58a) with b ββ , in order to obtain

where we used the following identity in the diffusive term:

Applying the superficial averaging operator (see Eq. (A1a) in Appendix A) corresponding to the β-phase to Eqs. (62) leads to

The last term on the above equation can be readily dropped because both the porosity and c ββ can be taken out from the averaging operator, leaving only b ββ . Therefore, on the basis of Eq. (58f), this last term is zero. Directing the attention to the first term on the rhs of the above equation, and applying the spatial averaging theorem, yields

As already noticed in this work, since closure variables are periodic, there are no spatial variations of average quantities.

After application of the interfacial boundary conditions given in Eq. (58d), this last result can be consequently expressed as follows:

The first term on the rhs of this equation is present in the definition of K ββ as indicated in Eq. (57a). Note that in the second term, the boundary condition given in Eq. (58c) was used. Moreover, since this term involves the closure variable b σβ , it is opportune to make the outer product of Eq. (58b) with b σβ and then apply the superficial averaging operator corresponding to the σ-phase, in order to obtain

Following developments similar to those used in the β-phase, it is not hard to deduce that

Using this expression in Eq. (66) and substituting the result into Eq. (64) lead to

APPENDIX B: TENSOR IDENTITIES

In this manuscript, many identities involving tensor algebra are used that are listed in the present appendix. For the sake of clarity, basic definitions and necessary material to recover them are first recalled in Gibbs' notation together with the Einstein summation convention. Let A and b be respectively a second-order tensor and a vector having all the necessary properties of regularity in the three dimensions of space for the derivation operators to be defined on their fields; these basic definitions are given by (we use the notation

Note here that ∇A is a third-order tensor,

Throughout the article, the nested convention for inner products is employed. This means that the indices that are closest together are those on which summation applies. As a consequence, the double inner product between two secondorder tensors A and B is defined as

Similarly, if C is a third-order tensor, we have

and further, if F is another third-order tensor

Importantly, while dealing with third-order tensors, one must be clear about the transposes of such algebraic quantities as three different transpose operators can be defined. In this manuscript, two of them are employed. The transpose denoted by the superscript T 1 permutes the first and second indices, whereas the transpose denoted by the superscript T 3 permutes the first and third indices, namely (C is again a third-order tensor),

In addition, throughout this work, the outer product between two vectors, a and b, is defined as

It is worth mentioning that some authors 52 denote this product as a ⊗ b. The outer product between a second-order tensor A and a vector b is given in a similar manner by

With this at hand, it is not hard to deduce the following formulas that are used in the article and which can be listed as follows.

(1) For a second-order tensor, A, which is a regular enough function of position, the following identities are applicable (the superscript T represents the classical transpose of a second-order tensor, i.e., A T ij = A ji ),

(2) For a second-order tensor A and a vector b, which are both continuous and derivable functions of position, the following identities are applicable:

When b is a solenoidal vector field, the following identities are valid:

(3) For any arbitrary second-order tensors A and B having the required regularity properties, the following identities hold:

∇A : ∇B T = ∇B T 3 : ∇A T 3 , (B23)

The last part of this appendix is dedicated to the proof that, for any second-order tensor A having the required regularity properties, the second-order tensor given by the expression B = ∇A T 3 : ∇A + ∇A T 1 is symmetric. The ij component of this tensor is given by

The ji component of this tensor is obtained by inverting the i and j indices, yielding

The first term on the rhs of both Eqs. (B25) and (B26) is obviously the same. Moreover, for the second term, one can notice that the two indices on which summation is performed are dummy and can hence be interchanged so that the second term on the rhs of both the above expressions is also the same, completing the proof that ∇A T 3 : ∇A + ∇A T 1 is a symmetric tensor.