Didier Lasseux
email: didier.lasseux@u-bordeaux.fr

Francisco J Valdés-Parada

On the developments of Darcy's law to include inertial and slip effects

Keywords:

The empirical Darcy law describing flow in porous media, whose convincing theoretical justification was proposed almost 130 years after its original publication in 1856, has however been extended to account for particular flow conditions. This article reviews historical developments aimed at including inertial and slip effects (respectively, when the Reynolds and Knudsen numbers are not exceedingly small compared to unity). Despite the early empirical extensions to include inertia and slip effects, it is striking to observe that clear formal derivations of physical models to account for these effects were reported only recently.

Introduction

Flows in porous media are of interest in numerous applications ranging from hydrology, hydrocarbon recovery, gas and nuclear waste storage, to drying of wood, transfer in food products or in living tissues to cite but a few. The main characteristic of this particular domain of fluid mechanics lies in the (sometimes extreme) complexity of the geometry of the channels where the flow takes place. Additionally, in many situations, this geometry is unknown in its very details and may vary over more or less long distances characteristics of heterogeneities. Within this context, the physical description of the flow in such materials may appear to be a tremendous challenge. 1 This certainly explains why empiricism remained so strong and lasted longer than in many other fields of fluid mechanics. In many situations, the interest is not in the details of the flow within the pores but rather in the flux-to-force governing laws at length scales including a large numbers of pores, although comprehensive analyses at the pore scale remain the corner stone in any progress towards the derivation of governing laws at larger scales. Clearly, active research in the description of transfer in porous materials was triggered by the publication of Darcy's law in the middle of the 19th century and the emergence of a key macroscopic physical characteristic of a porous medium, namely the ability of a fluid to flow through it, i.e. its permeability.

Darcy's law as an empirical basis

Ever since its empirical formulation in 1856, Darcy's law [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon[END_REF] has been a hallmark in modeling momentum transport through porous media. In this classical publication, there is a section entitled Détermination des lois d'écoulement de l'eau à travers le sable, dedicated to the study of water flows through a bed of sand where the following relation is proposed (see page 594 in [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon[END_REF]): q = k s e (h + e) [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon[END_REF] where q is the volumetric flow rate, s is the cross section of the sand bed, e is the bed width, h is the pressure (or head) difference between the surface and the base of the sand bed, and finally k was proposed as a coefficient that depends on the permeability of the bed 2 and on the properties of the fluid. For an excellent review about the origin of Darcy's law, the interested reader is referred to the work by Zerner [START_REF] Zerner | Aux origines de la loi de darcy[END_REF]. The use of this simple relation requires that the only resistance to the flow through the porous medium is due to viscous stresses induced by an isothermal, creeping (or laminar) steady flow of a Newtonian fluid within an inert, rigid and homogeneous porous medium. However, the lack of a rigorous upscaling technique prevented a formal derivation of this equation until the late 1960s, as it will be detailed later.

For a very long period of time -around fifty years -this law has been essentially verified experimentally in its global form, but was not considered in a local differential form nor derived on a theoretical basis. One finds a differential expression in the analysis of the flows in aquifers by J. Boussinesq [START_REF] Boussinesq | Recherches théoriques sur l'écoulement des nappes d'eau infiltrées dans le sol et sur le débit des sources[END_REF] as a result of an analogy with heat transfer in a continuum. This work also reports an extension of the flow-rate-to-head-gradient relationship to non-homogeneous media. A formal derivation of a 1D local expression of this law obtained from the solution to the Stokes equation for a flow parallel to a regular array of infinite parallel cylinders (sufficiently apart from each other, i.e. for relatively large porosities) is due to Emersleben in 1925 [START_REF] Emersleben | The Darcy filter formula[END_REF]. A derivation mainly based on dimensional analysis was later proposed by Muskat and Botset in 1931 [START_REF] Muskat | Flow of gas through porous materials[END_REF] for a compressible flow in which the pressure difference is recognized to be replaced by the difference of the squares of the pressures. Substantial literature will then appear during the 1950s, in which many different approaches to demonstrate Darcy's law will be tested (see for instance [START_REF] Irmay | On the theoretical derivation of Darcy and Forchheimer formulas[END_REF] and references therein). Although these articles helped progressing into the understanding of the applicability of Darcy's law, almost all of them relied on analogies, hypotheses or postulates that left them incomplete. The first extension to three-dimensions and to non-isotropic materials was reported by Hall [START_REF] Hall | An analytical derivation of the Darcy equation[END_REF], who introduced a permeability tensor, which is also based on some pre-requisites (see in particular Eq. ( 17) therein and the way the permeability is identified).

Despite the lack of formal derivation of Darcy's law, which can be expressed for a 1D flow in the x-direction as [START_REF] Muskat | Flow of Homogeneous Fluids through Porous Media[END_REF] q = -

K s μ ∂ p β β ∂x ( 2
)
the meaning of the permeability and its relationship to the underlying pore structure focused closed attention. In the above expression, K is the permeability having units of m 2 and μ the fluid viscosity. In addition, p β β is the intrinsically-averaged pressure in the porous medium, defined as:

p β β = 1 V β V β p β dV (3)
Here, V β is the domain (of volume V β ) occupied by the fluid phase β within a representative averaging domain (or REV) (see Fig. 1), and p β is the pore-scale pressure. An early estimate of K was inspired by an analysis due to Blake in 1922 [START_REF] Blake | The resistance of packing to fluid flow[END_REF] of flow over packings of grains of different shapes and a comparison to flows in capillary tubes that resulted in the following estimate

K = 1 k 0 S 2 0 ε 3 (1 -ε) 2 (4)
where ε is the porous medium porosity, S 0 denotes the specific surface of the particles and it is defined as the ratio of the area of the particle to its volume. The coefficient S 0 , related to the effective particle diameter, d p , was identified from an analogy with spherical particles by

S 0 = 6 d p ( 5
)
2 In [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon[END_REF], H. Darcy indicated that k "depends on the permeability of the sand layer". In fact, k is the hydraulic conductivity, having the unit of a length per Finally, k 0 is an adjustable coefficient that was later known as the Kozeny [START_REF] Kozeny | Ueber kapillare leitung des wassers im boden[END_REF] coefficient. Carman [START_REF] Carman | Flow of Gases Through Porous Media[END_REF] suggested taking k 0 = 5, so that equation ( 4) can now be expressed as

K = d 2 p ε 3 180(1 -ε) 2 (6)
which is usually known as the Kozeny-Carman equation and which is often used, sometimes abusively, to predict the value of permeability. Certainly, a precise correlation for any type of pore structure is out of reach.

Theoretical foundation of Darcy's law

The lack of convincing formal derivations of Darcy's law that lasted over more than a century is obviously related to the lack of clear upscaling methods allowing one to obtain the macroscopic conservation equations form their microscopic (i.e. pore-scale) counterparts, as pointed out by Zerner [START_REF] Zerner | Aux origines de la loi de darcy[END_REF]. These methods emerged in the 1970s and a first attempt was proposed by S. Whitaker in 1966 [START_REF] Whitaker | The equations of motion in porous media[END_REF], who clearly obtained a generalization of Darcy's law in the following vectorial form: [START_REF] Hall | An analytical derivation of the Darcy equation[END_REF] where ρ is the fluid density, g is the gravity acceleration, K is the permeability tensor and v β the seepage velocity, which is defined as the superficial average:

v β = - K μ • ∇ p β β -ρg
v β = 1 V V β v β dV (8)
with V being the volume of the REV and v β the pore-scale velocity vector. The same result was achieved exactly in the same period by C. Marle [START_REF] Marle | Écoulements monophasiques en milieu poreux[END_REF]. However, in these articles, no clear structural link (i.e. a closure) is provided between the micro-and the macroscale allowing one to infer the dependence of the permeability upon the pore-scale geometry. It was not earlier than in the middle of the 1980s that a closed rigorous form was achieved by Whitaker [START_REF] Whitaker | Flow in porous media I: a theoretical derivation of Darcy's law[END_REF] using the volume averaging method [START_REF] Whitaker | The Method of Volume Averaging[END_REF], which included an intrinsic closure scheme allowing one to predict the values of the components of the permeability tensor K. Some derivations of Darcy's law have used other upscaling techniques such as homogenization [START_REF] Auriault | Nonsaturated deformable porous media: quasistatics[END_REF]. These and other techniques depart from the governing equations at the pore scale and, after the application of an averaging operator (such as the one sketched in Fig. 1) and many mathematical steps, lead to an upscaled model in terms of effective-medium coefficients that capture the essential (i.e. non-redundant) information from the pore scale. In this way, the permeability tensor can be viewed as a signature of the porous medium topology at a scale that is larger than the pore scale.

Over the past century, there have been some modifications brought to Darcy's law that extended its applicability to more complicated transport processes than those originally considered by Darcy. Among the extensions to Darcy's law, a non-exhaustive list should include: 1) Forchheimer's modification to allow for the study of non-creeping flow regimes [START_REF] Forchheimer | Wasserbewegung durch boden[END_REF]; 2) Brinkman's correction to include macroscopic viscous stress by introducing an effective viscosity coefficient [START_REF] Brinkman | A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles[END_REF]; 3) Klinkenberg's modification of the permeability tensor to study gas slip flows in porous media [START_REF] Klinkenberg | The Permeability of Porous Media to Liquids and Gases[END_REF]; 4) application to heterogeneous media by means of large-scale volume averaging [START_REF] Quintard | Écoulement monophasique en milieu poreux: effet des hétérogénéités locales[END_REF]; 5) non-isothermal flow of non-Newtonian fluids in porous media (cf. [START_REF] Kairi | Effect of viscous dissipation on natural convection heat and mass transfer from vertical cone in a non-Newtonian fluid saturated non-Darcy porous medium[END_REF] for instance).

In what follows, the focus is laid upon two extensions to Darcy's law that are of major importance, namely the inertial one-phase flow and the gas slip flow in homogeneous porous media. Our aim is to carefully review these extensions and draw some conclusions and perspectives on these subjects.

Inertial one-phase flow in porous media

The analogy with flows in ducts has been widely used in the derivation of empirical flow models in porous media and might have been a source of inspiration for H. Darcy to obtain the filtration law he reported in his book [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon[END_REF]. As mentioned in a remarkable analysis by Zerner [START_REF] Zerner | Aux origines de la loi de darcy[END_REF], H. Darcy dedicated specific experiments to verify Poiseuille's law in the context of slow flow. One of the major motivations was his questioning of the validity of the relationship between the "pressure drop" P and the average velocity u in a tube of length L, which was then of common use, i.e.

P L = au + bu 2 (9)
a relationship mainly due to du Buat and Gaspard Riche de Prony (see [START_REF] Zerner | Aux origines de la loi de darcy[END_REF]), with a and b being coefficients that had to be determined experimentally.

Forchheimer's correction to Darcy's law

It is striking to observe that the form of Eq. ( 9) corresponds to the relationship proposed by Forchheimer [START_REF] Forchheimer | Wasserbewegung durch boden[END_REF] to account for "rapid" flows in porous media with the classical analogy on u taken as the filtration or seepage velocity in 1D:

P L = μ K u + ρβu 2 (10)
where ρ is the fluid density and β is the Forchheimer coefficient, also called the coefficient of inertial resistance or inertial resistance factor.

The quadratic correction introduced by Forchheimer about half a century after the empirical validation and publication of Darcy's law was obviously inspired from Eq. ( 9), despite this remarkable time lapse between Forchheimer's publication and that of Darcy. In 1931, Muskat and Botset [START_REF] Muskat | Flow of gas through porous materials[END_REF] reported experimental results of gas flow through different types of porous materials, in which they observed that the gradient of the square of the pressure was proportional to a power of the mass velocity ranging between 1 and 2 (1 when the flow was "completely viscous", and 2 when it was "completely turbulent"). This empirical form of Eqs. ( 9) and ( 10) was accepted to account for inertial flows in porous media over an additional half-century during which only few alternative forms were put forth, like for instance [START_REF] Polubarinova-Kochina | Theory of Ground Water Motion[END_REF] P L = au + bu 2 + cu 3 [START_REF] Carman | Flow of Gases Through Porous Media[END_REF] with a, b and c being adjustable coefficients. During the fifty years following Forchheimer's publication, some confusion remained about the physical origin of the quadratic correction to Darcy's law as it was often attributed to turbulence, although some references made clear statements on that point [START_REF] Irmay | On the theoretical derivation of Darcy and Forchheimer formulas[END_REF]. In fact, Irmay [START_REF] Irmay | On the theoretical derivation of Darcy and Forchheimer formulas[END_REF] argued that there is no reason, in general, to expect a linear solution to the non-linear Navier-Stokes equations. The Poiseuille solution in straight tubes is an exception caused by vanishing curvature, which is not the case in real porous media. Agreement has been quite unanimous on the threshold value of the Reynolds number at which the correction becomes significant, i.e. for 1 ≤ Re d ≤ 15, where

Re d = ρ v β d g μ
is based on the filtration velocity v β and the typical grain size d g of the porous material. However, it was not before 1962 and the publication by Chauveteau and Thirriot [START_REF] Chauveteau | Régimes d'écoulement en milieu poreux et limite de la loi de Darcy[END_REF], in which a flow regime classification was proposed, that turbulence in this range of Reynolds numbers was dismissed. Turbulence has been confirmed to typically arise for Re d ∼ 100 [START_REF] Dybbs | A new look at porous media fluid mechanics -Darcy to turbulent[END_REF][START_REF] Ghaddar | On the permeability of unidirectional fibrous media: a parallel computational approach[END_REF][START_REF] Koch | Moderate Reynolds number flows through periodic and random arrays of aligned cylinders[END_REF].

During this period, and even up to the end of the 1970s, comparisons with experimental results were reported, showing quite good agreement for various types of porous structures, including packed beds of grains, bundles of capillary tubes or fibrous media for flows of gases or liquids [START_REF] Ergun | Fluid flow through packed columns[END_REF][START_REF] Dullien | Flow rate-pressure gradient measurement in periodically nonuniform capillary tubes[END_REF]. From a practical point of view, the Forchheimer model has been used for applications in hydrology, petroleum and chemical engineering [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF][START_REF] Geertsma | Estimating the coefficient of inertial resistance in fluid flow through porous media[END_REF]. From a theoretical point of view, the same period was marked by various attempts to justify the form of Eq. ( 9). This was carried out on the basis of different approaches ranging from simplified derivations from the Navier-Stokes equations or analogies with pressure drop through capillary orifices [START_REF] Irmay | On the theoretical derivation of Darcy and Forchheimer formulas[END_REF]. The emergence of more formal upscaling techniques during the 1970s has led to further developments in the following years that also attempted to justify the quadratic form of the inertial correction to Darcy's law [START_REF] Marle | Écoulements monophasiques en milieu poreux[END_REF][START_REF] Hassanizadeh | High velocity flow in porous media[END_REF][START_REF] Giorgi | Derivation of the Forchheimer law via matched asymptotic expansions[END_REF][START_REF] Chen | Derivation of the Forchheimer law via homogenization[END_REF].

Refinements on the inertial correction to Darcy's law

Impulsed by the development of both computational resources and numerical methods, the velocity dependence of the corrective term on Darcy's law regained much attention from the early 1990s on. It was during this period that the quadratic velocity dependence of this correction was questioned and analyzed in depth. Numerical simulations carried out h This model, subject to time (μt * ρ 2 β ), length-scales ( β r 0 L) and Reynolds number

ρ v β β μ β L
through networks of parallel cylinders of circular cross-section for a flow orthogonal to the cylinders showed that the correction scales as a 3rd power of the filtration velocity [START_REF] Barrère | Modélisation des écoulement de Stokes et de Navier-Stokes en milieu poreux[END_REF] instead of a quadratic term. This behavior was theoretically confirmed quasi simultaneously, regardless the geometry at the pore-scale, from formal derivations based on a rigorous upscaling procedure using double-scale homogenization with a closure process involving periodic representations of the porous medium [START_REF] Wodie | Correction non linéaire de la loi de Darcy[END_REF][START_REF] Mei | The effect of weak inertia on flow through a porous medium[END_REF][START_REF] Rasoloarijaona | Nonlinear seepage flow through a rigid porous medium[END_REF]. The exponent 3was confirmed for a Reynolds number, Re p , based on the characteristic pore size ranging between δ 1/2 and 1, where δ is the micro-to-macroscale ratio. This result was further emphasized later on [START_REF] Skjetne | High velocity flow in a rough fracture[END_REF] and, from its original evidence, led to identify this regime as the so-called "weak inertia regime", formalized for homogeneous isotropic and periodic media. Additional numerical works over larger intervals of the Reynolds number for flows in many different structures confirmed the existence of weak inertia and extended the classification of flows under (at least) three distinct regimes with crossovers, namely [START_REF] Skjetne | New insights on steady, non-linear flow in porous media[END_REF] i) the weak inertia regime occurring at the onset of non-linearity in the flow-rate-to-pressure-drop relationship for δ 1/2 ≤ Re p ≤ 1; ii) the "strong inertia regime" characterized by a correction to Darcy's law that scales as the square of the filtration velocity, i.e. leading to a Forchheimer type of model for Reynolds numbers in the range 1to 10; iii) the turbulent regime appearing for Reynolds numbers typically of the order of 100. Nevertheless, a detailed physical explanation of a cubic or quadratic correction to Darcy's law in order to account for inertia still leaves much to be desired and remains a widely open question [START_REF] Skjetne | High-velocity laminar and turbulent flow in porous media[END_REF], while the description of the non-linearity is essentially qualitative. Obviously, inertial macroscopic forces cannot be invoked as they remain negligible compared to viscous forces [START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF], and this can be proved to hold as long as Re p δ -1 . Consequently, the nature of the non-linearity in the relationship between the macroscopic drag force and filtration velocity must certainly be explained from the signature of viscous and inertial forces at the microscale. On the one hand, several mechanisms suggesting that inertia alone, at the pore-scale, can explain the macroscale behavior may be put forth such as: i) streamlines bending due to the tortuosity of the structure and to local converging-diverging flow patterns; ii) backflows and separations resulting from form drag; iii) pore networks actively involved in the flow that are velocity-dependent as a consequence of i) and ii) yielding variations in the dissipation of the kinetic energy [START_REF] Andrade | Inertial effects on fluid flow through disordered porous media[END_REF]. On the other hand, microscale viscous drag effects may be considered to contribute to the non-linearity when boundary layers at the solid-fluid interfaces, which become thicker when the Reynolds number increases (see an experimental observation in [START_REF] Dybbs | A new look at porous media fluid mechanics -Darcy to turbulent[END_REF]), are taken into account. With this mechanism, the inertial core flow (in the center of the pores) may be easily understood as being strongly dependent upon the Reynolds number, partly explaining the different regimes.

Further developments

Although the existence of the two regimes (weak and strong inertia) has been widely accepted, a lot still needs to be understood regarding the universal existence and dependence of these regimes (and the associated crossovers) upon many parameters such as porosity, structural order, anisotropy, etc. In addition, most of the experimental or numerical characterizations of the macroscale inertial correction have been carried out in 1D in a scalar form although many references pointed out that tensorial coefficients must be involved [START_REF] Hassanizadeh | High velocity flow in porous media[END_REF][START_REF] Mei | The effect of weak inertia on flow through a porous medium[END_REF][START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF][START_REF] Giorgi | Derivation of the Forchheimer law via matched asymptotic expansions[END_REF]. Even if the development did not allow for the formal identification of the above-mentioned regimes, a convincing derivation, relying on rigorous upscaling, of a macroscopic model for momentum transport of one-phase flowing in homogeneous porous media with inertia is certainly due to W itaker [START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF].

1 constraints that were clearly formulated, provides the macroscale momentum equation, which reads [START_REF] Whitaker | The equations of motion in porous media[END_REF] or equivalently

v β = - H μ • ∇ p β β -ρg
v β = - K μ • ∇ p β β -ρg -F • v β (13)
where v β is the filtration (or seepage) velocity as defined in Eq. ( 8), H is the apparent permeability, K the intrinsic permeability and F the inertial correction tensor (which is a function of the filtration velocity); ∇ p β β is the macroscopic pressure gradient, p β β being defined in Eq. (3). Details on the averaging method employed to obtain this result are given in [START_REF] Whitaker | The Method of Volume Averaging[END_REF]. In addition to the macroscopic model, the upscaling provides the means to determine the associated macroscopic coefficients (K, H and F contained in Eqs. ( 12) and ( 13)) from the solution to the ancillary problems (so-called closure problems). Undoubtedly, while completing the physical description, this contribution concludes some previous derivations on the same problem performed with homogenization [START_REF] Mei | The effect of weak inertia on flow through a porous medium[END_REF][START_REF] Skjetne | High-velocity laminar and turbulent flow in porous media[END_REF] and opens new perspectives to more in-depth investigation of the inertial correction to Darcy's law. In this spirit, closure problems were solved in order to compute H and F for different ordered and disordered model structures so as to investigate the existence of the regimes and their dependence upon the porosity and the microscale pore structure [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF]. The main important conclusions deriving from this analysis can be summarized as follows. The two tensors H and F are generally dense and non-symmetric for ordered structures, even if the medium is isotropic at the macroscale in the Darcy regime. The non-symmetry coincides with a macroscopic drag force, which is not necessarily aligned with the macroscopic velocity. Symmetry is recovered for specific pressure gradient orientations along symmetry axes of the structure (when present). Dissymmetry of these tensors decreases when structural disorder increases. Whatever the structure under concern, the weak inertia regime (a cubic inertial correction to Darcy's law) is always observed. For ordered structures, the strong inertia regime (i.e. aquadratic (Forchheimer) correction to Darcy's law) does not necessarily exist and it is otherwise restricted to a narrow interval of the Reynolds number. For disordered structures, the Forchheimer type of correction is a robust approximation over a very significant range of values of the Reynolds number. Moreover, the crossover value of the Reynolds number at which the quadratic correction becomes relevant decreases with increasing structural disorder, the weak inertia regime being restricted to a small range of Reynolds numbers where the correction is not very significant. This certainly explains why this regime is most of the time overlooked in experimental investigations on porous media having a random pore structure. In any situation, using a model such as that proposed in Eq. ( 9) implies that the permeability in the linear term in the filtration velocity differs from the intrinsic permeability. Even if some progress has been achieved in the physical explanation and the theoretical derivation of formal models to account for inertia effects for one-phase flows in porous media, much work remains to be done to understand fully the non-linearities associated with this type of process. For instance, most of the analyses so far were dedicated to the laminar steady regime and the occurrence of unsteadiness remains widely unexplored (this problem was barely outlined recently in [START_REF] Agnaou | From steady to unsteady laminar flow in model porous structures: an investigation of the first Hopf bifurcation[END_REF]), as well as the turbulent regime; however, further discussion of these topics is out of the scope of the present review.

Slip flow in porous media

Gas flows in porous media differ considerably from liquid-phase flows, in particular for situations in which the pore sizes are comparable to the mean free path of the gas molecules. This is the case in many practical applications including microand nano-fluidic systems such as MEMS and nano-porous media, transport in fibrous media, gas flow during soil remediation, long-term nuclear waste disposal, among many others. Due to the current high relevance of this type of transport, the evolution from the pioneering works performed in the nineteenth century to some of the current developments are briefly summarized in this section.

Slip flow background

A one-phase flow in confined systems of dimensions comparable to the mean free path leads to rarefaction effects that give rise to many interesting contributions in transport phenomena. In his classical study of the stresses that rarefied gases experience, J.C. Maxwell [START_REF] Maxwell | On stresses in rarified gases arising from inequalities of temperature[END_REF] proposed that, close to the surface of a solid, there should be a sliding of the gas in contact with the solid in the direction of a tangential stress. Maxwell proposed that the velocity should be proportional to the tangential stress and inversely proportional to the viscosity of the fluid. Under isothermal conditions, the sliding velocity for a 1D flow in the x-direction is expressed as:

v = G dv dx ( 14
)
where the coefficient G is the coefficient of slipping, defined as

G = 2 3 2 f -1 λ β ( 15
)
with λ β being the mean free path and f the fraction of gas molecules that are diffusely scattered at the surface. Hence, if the solid surface is wholly absorbent, G = 2/3λ β . The proposal from Maxwell is consistent with a previous study by Navier [START_REF] Navier | Mémoire sur les lois du mouvemet des fluides[END_REF],

and we shall thus refer to Eq. ( 14) as the Navier-Maxwell equation. It should be noted that Maxwell obtained Eq. ( 14) considering a single-component gas; the extension of this equation to multicomponent mixtures is found in Jackson [START_REF] Jackson | Transport in Porous Catalysts[END_REF].

During the last quarter of the nineteenth century, rarefaction effects were shown to increase significantly the flow rate with respect to that predicted from Poiseuille's law [START_REF] Kundt | On friction and heat-conduction in rarefied gases[END_REF]. This motivated several investigations during the early years of the twentieth century, in particular those by M. Knudsen. The scattering of gas molecules from solid walls was fundamental in Knudsen's theory [START_REF] Knudsen | Die gesetze der molekularstromung und der inneren reibungsströmung der gase durch röhren[END_REF]. Knudsen studied molecular flows in tubes and determined the dependence on tube dimensions. He discussed the transition from Poiseuille's flow in terms of the ratio of the mean free path of the gas molecules to the characteristic size of the apparatus (say, β ). The reason for considering this ratio is due to the fact that, in the continuum approach, the slip velocity may be understood as the average flow velocity of the molecules at a distance from the wall that is equal to the mean free path [START_REF] Dullien | Porous Media. Fluid Transport and Pore Structure[END_REF]. Hence, as the mean free path becomes a bigger fraction of the tube diameter, the slip velocity increases with respect to the bulk velocity. This important ratio between the mean free path and the tube diameter is nowadays known as the Knudsen number in his honor (i.e. K n = λ β / β ). In the comprehensive review by Steckelmacher [START_REF] Steckelmacher | Knudsen flow 75 years on: the current state of the art for flow of rarefied gases in tubes and systems[END_REF], the historical development and relevance of Knudsen's works are presented in detail.

Later on, Adzumi [START_REF] Adzumi | Studies on the flow of gaseous mixtures through capillaries. I. The viscosity of binary gaseous mixtures[END_REF][START_REF] Adzumi | Studies on the flow of gaseous mixtures through capillaries. II. The molecular flow of gaseous mixtures[END_REF][START_REF] Adzumi | Studies on the flow of gaseous mixtures through capillaries. III. The flow of gaseous mixtures at medium pressures[END_REF], published a series of papers dedicated to the study of gas flow through capillaries. He considered three cases: 1) when λ β is very small in comparison to the diameter of the capillary (i.e. K n 1); 2) when λ β is large compared to the diameter (i.e. K n > 1), and 3) when λ β is comparable to the capillary diameter (i.e. K n ∼ 1). In the first case, he found Poiseuille's law to be quite suitable, and the flow is inversely proportional to the gas viscosity [START_REF] Adzumi | Studies on the flow of gaseous mixtures through capillaries. I. The viscosity of binary gaseous mixtures[END_REF]. In the second case, the flow rate was found to be independent of the viscosity, but inversely proportional to the square root of the gas molecular weight, M, [START_REF] Adzumi | Studies on the flow of gaseous mixtures through capillaries. II. The molecular flow of gaseous mixtures[END_REF]. The flow characteristics in the third case were found by Adzumi to be a combination of the two first ones [START_REF] Adzumi | Studies on the flow of gaseous mixtures through capillaries. III. The flow of gaseous mixtures at medium pressures[END_REF]. Nowadays, there is some consensus that the following bounds are identifiable: for Kn < 10 -3 , the laws of continuum mechanics are safely applicable, and non-slip can be assumed at the solid boundaries; for 10 -3 < Kn < 10 -1 , the flow regime corresponds to slip flow, and the Navier-Maxwell equation must be considered at the solid-fluid interface;

for 10 -1 < Kn < 10, there is a transition regime in which the laws of continuum mechanics are likely to fail because the continuum hypothesis is no longer satisfied. For porous media applications in the transition regime, Maxwell [START_REF] Maxwell | II. Illustrations of the dynamical theory of gases[END_REF] proposed that the action of the porous material over the gas was similar to a number of dust particles of the moving system, hence giving rise to the well-known dusty gas model. This model has the nice feature that the interactions of gas molecules with the dust molecules simulate their interaction with the rigid porous matrix, thus avoiding the problem of flux variations across the sections of the pores [START_REF] Jackson | Transport in Porous Catalysts[END_REF]. Finally, for Kn > 10 the gas kinetic theory must be considered because description in terms of particle-wall collision operators is required. This limiting situation is called molecular streaming or Knudsen flow, and it is characterized by the fact that the flow takes place by diffusion, instead of viscous, mechanisms [START_REF] Dullien | Porous Media. Fluid Transport and Pore Structure[END_REF]. Here, the term diffusion means that the flow results from creep at the wall rather than from molecule-to-molecule collisions.

Pore-scale slip flow and its consequence on Darcy's law

Few years after Adzumi's works, in his study of the permeability of gas flows in porous media in the slip regime, Klinkenberg [START_REF] Klinkenberg | The Permeability of Porous Media to Liquids and Gases[END_REF] found that this coefficient is almost a linear function of the reciprocal mean pressure, p β β . Using an idealized porous medium representation consisting in an array of capillaries, Klinkenberg was able to deduce the following relation between the apparent permeability, K s , and the intrinsic permeability, K :

K s = K 1 + b p β β ( 16
)
with b being a surface-and gas-dependent constant. Evidently, at sufficiently large gas pressures, K s approaches K . Klinkenberg used the above equation to predict the values of K s /K in different experimental conditions with reliable accuracy, thus

showing the need to consider slip effects in the determination of the permeability. In this way, the Darcy-Klinkenberg model is (gravity is omitted here):

v β = - K μ 1 + b p β β ∂ p β β ∂x (17)
for an average one-dimensional flow in the x-direction, while μ is the fluid dynamic viscosity.

An important point of discussion in the recent literature related to gas transport in porous media is about the pertinence of the linear and first-order interfacial boundary condition given in Eq. ( 14). Shen et al. [START_REF] Shen | A kinetic-theory based first order slip boundary condition for gas flow[END_REF] derived a first-order slip condition from the Chapman-Enskog solution to the Boltzmann equation. This approach consists in linearizing the Boltzmann equation using a perturbation expansion for the probability function in terms of the Knudsen number. The resulting condition includes an additional term due to the pressure gradient along the flow's direction. The success of first-order models has been argued to be constrained to slightly rarefied gas flows by Deissler [START_REF] Deissler | An analysis of second-order slip flow and temperature-jump boundary conditions for rarefied gases[END_REF]. According to this author, as the pressure in the gas becomes smaller, the velocity profiles may be nonlinear over a distance from the solid surface corresponding to the mean free path and the jumps at the interface may be expected to be functions of higher-order normal and tangential derivatives. Deissler thus proposed to use a second-order boundary condition that matches the Navier-Stokes equations for slip flows, finding good agreement with experimental results. Another early proposal of second-order models is the scalar one by Cercignani [START_REF] Cercignani | Higher Order Slip According to the Linearized Boltzmann Equation[END_REF] on the basis of the Bhatnagar-Gross-Krook (BGK) approximation of the Boltzmann equation. However, these extensions consist of modifications to the scalar equation [START_REF] Whitaker | Flow in porous media I: a theoretical derivation of Darcy's law[END_REF], in other words, they are constrained to simple geometries where one-dimensional flow is applicable. Unfortunately, extensions to more complicated geometries are not straightforward and they remain a challenge. Furthermore, since the Navier-Stokes equations are first-order accurate in the Knudsen number, it is not easy to justify the use of higher-order boundary conditions. Hence, an alternative is to use higher order momentum transport equations, such as the Burnett equation [START_REF] García-Colín | Beyond the Navier-Stokes equations: Burnett hydrodynamics[END_REF], which are the result of keeping the second-order terms in the Chapman-Enskog approach to approximate the Boltzmann equation.

As mentioned above, one of the limitations of the Navier-Maxwell boundary condition is that it is restricted to onedimensional flows. This limitation is also shared by the Klinkenberg model. To address this issue, Einzel et al. [START_REF] Einzel | Boundary condition for fluid flow: curved or rough surface[END_REF] proposed a generalized version of the slip boundary condition, which can be expressed as follows:

v β = - 2 -σ ν σ ν ξ λ β n • ∇v β + ∇v T β • (I -nn) (18)
where σ ν is the tangential-momentum accommodation coefficient, n is the unit normal vector directed from the fluid toward the sold phase and I is the identity tensor. The coefficient σ ν accounts for the average tangential momentum exchange between the molecules and the fluid, and can vary from zero to one.

Theoretical macroscopic slip flow models in porous media

In the same line of thoughts as those mentioned for inertial flows, the Klinkenberg model, although widely used in the literature, was formally derived only more than half a century after its publication. Using the homogenization method, Skjetne and Auriault [START_REF] Skjetne | Homogenization of wall-slip gas flow through porous media[END_REF] considered the steady-state, low-velocity Navier-Stokes equations for compressible flows in the slip regime, i.e. Re Kn 1(for Re = O(δ = /L), with and L being the characteristic length scales at the microscale and macroscale, respectively). In this work, for the first time, a vectorial form of the Klinkenberg model was rigorously deduced for conditions in which both local compressibility and inertial effects are smaller than the wall-slip effects. The apparent permeability tensor was found to be positive-definite and non-symmetric, in general. This study was subsequently expanded by Chastanet et al. [START_REF] Chastanet | Flow of low pressure gas through dual-porosity media[END_REF] to derive the corresponding upscaled model for low-pressure gas flows in dual-porosity media including fractures. Their study highlighted that the Knudsen number should be considered in addition to the separation of characteristic length scales in the system in order to assess the domains of validity of upscaled models for gas flow.

The derivation of the effective-medium equation corresponding to slightly-compressible slip-flow conditions using the volume averaging method has been carried out recently by Lasseux et al. [START_REF] Lasseux | A macroscopic model for slightly compressible gas slip-flow in homogeneous porous media[END_REF]. This work completes the previous derivations by Skjetne and Auriault on the following points: 1) the compressibility effects are taken into account in the framework of slightly compressible flows restricted to small Reynolds numbers and small frequency number; 2) the vectorial form of the slip boundary condition is considered in its complete form, i.e. including the complete shear-rate at the fluid-solid interface as shown in Eq. ( 18); 3) it is derived for a barotropic fluid without any assumption on the equation of state of the gas. For an ideal gas, the resulting upscaled model is

v β = - 1 μ K • ⎛ ⎝ I + ξ μ p β β πR T β β 2M S ⎞ ⎠ (19)
which involves two tensors, namely, the intrinsic permeability tensor, K, and a slip-flow correction tensor, S. The parallelism between Eqs. ( 17) and ( 19) is obvious. The ancillary closure problem required to predict the values of the effective-medium coefficients was derived and formally solved for simple porous medium geometries in two-and three-dimensional unit cells. Furthermore, the dependence of s (S = sI) on K (K = K I) was found to obey a power-law relationship, with the value of the exponent depending on the geometrical configuration. In a subsequent work by Lasseux et al. [START_REF] Lasseux | An improved macroscale model for gas slip flow in porous media[END_REF], the slip correction was more accurately described by considering an expansion in the Knudsen number at the closure level. This leads to a reformulation of the closure problem as a differential (instead of an integro-differential) boundary-value problem, which was solved in more complicated unit cells in order to predict the apparent permeability tensor. Furthermore, with this expansion, the slip-flow correction tensor was shown to be the sum of slip corrections at the successive orders of the Knudsen number. The consideration of the complete form of the boundary condition at the solid-fluid interface was found to be crucial for the prediction of the slip corrections at the different orders of Kn. Their analysis evidenced a nonlinear relationship between the apparent permeability and the Knudsen number. This relationship motivates further theoretical and experimental research on the subject, in particular for highly porous structures.

Conclusions

This work has been dedicated to the analysis of the evolution of two major modifications to Darcy's law, namely the inclusion of inertial and slip effects. The review carried out for both extensions suggests the following conclusions and prospects.

-The empirical introduction of a quadratic correction in terms of the filtration velocity to Darcy's law by Forchheimer in 1901 for 1D flows was certainly inspired by a model that was commonly used for flows in pipes prior to the publication of Darcy's law. This correction has been widely supported by empiricism for more than 90 years, a period after which a first theoretical 3D model was achieved by upscaling (homogenization), showing that the onset of deviation from Darcy's law due to inertia involves a correction that rather scales as a 3rd power of the filtration velocity in a weak-inertia regime. For larger Reynolds number values, a so-called strong inertia regime, where the Forchheimer correction is due to hold, was accepted. An upscaled 3D complete model, obtained by volume averaging five years later, was used to highlight the fact that, if the weak inertia regime always exists, the quadratic correction does not hold in some particular situations of pore-scale ordered structures and that, however, in the presence of disorder, the Forchheimer-type correction is a robust one. At this point, a quite different observation from that indicated below for slip flows may be pointed out regarding corrections made to Darcy's law in order to include inertial effects. Indeed, if the underlying physics at the pore scale does not rise any particular question, issues are mainly related to the understanding of the different flow regimes at the macroscale, the physical mechanisms that trigger the transition from one regime to another and the associated range of Reynolds numbers, together with the possible occurrence of unsteadiness. -The original identification of the slip effect for gas flows in porous media, when the Knudsen number is not exceedingly small compared to unity, by [START_REF] Klinkenberg | The Permeability of Porous Media to Liquids and Gases[END_REF], which led this author to propose a correction to the intrinsic permeability (in the 1D case) inversely proportional to the mean pressure, has been accepted with empirical justifications for almost 60 years. After this period, a first theoretical 3D model was derived, followed by more refined ones 15 years later using rigorous upscaling techniques. While these upscaling tools are now at hand for such theoretical developments, it clearly appears that the main issues lie in an appropriate pore-scale description of the physics in that case, and in particular, in the slip-boundary condition, which still requires some efforts to better capture, and possibly extend, its domain of application in terms of the Knudsen number. In the same spirit, despite some attempts, macroscale models for transitional and strongly rarefied flow regimes still require important efforts. In parallel, numerical and experimental investigations are necessary to highlight the understanding of the detailed physics.

As a prospect, one may finally conclude with the following. From the evolution of two modifications to Darcy's law investigated in the present work, it is important to remark that, after one-dimensional corrections were proposed, followed by extensive experimental analyses, there was a long time period before rigorous deductions were presented. This delay can be explained by the fact that the different necessary theoretical frameworks for performing upscaling from the pore-scale to the macroscale are relatively new (about 30-40 years old) and also by the fact that, in many situations, the one-dimensional versions remained relatively satisfactory. However, with recent applications directed to micro-and nanofluidic devices, among others, there is a need for more accurate macroscale models, that can be validated through comparison with reliable experimental data and direct numerical simulations. In this last issue, there is a strong tendency in current research towards understanding macroscale phenomena through imaging and direct simulations. This has been made possible by recent advances in computational capabilities. With this perspective in mind, it is relevant to pose the question of which direction future advances of flows in porous media will take. Probably, experimental and numerical studies would still evolve very significantly. However, strong efforts should certainly be dedicated to more sophisticated upscaling approaches applicable to more challenging situations than those studied in the current state of the art.

Fig. 1 .

 1 Fig. 1. Sketch of a porous medium including an averaging domain and the phases involved.

The one-phase slow flow is probably one of the simplest mechanism one can think about and there is a lot of other tremendously more complex physical processes in porous media of relevance from both scientific and industrial points of view, including multiphase flows, compressible flows, phase change, deformable porous media, reactions in multicomponent systems, etc.


Catherine Teresa Paez-Garcfa

Francisco J Valdes-Parada

Didier Lasseux
email: didier.lasseux@u-bordeaux.fr

Macroscopic momentum and mechanical energy equations for incompressible single-phase flow in porous media

Modeling flow in porous media is usually focused on the governing equations for mass and momentum transport, which yield the velocity and pressure at the pore or Darcy scales. However, in many applications, it is important to determine the work (or power) needed to induce flow in porous media, and this can be achieved when the mechanical energy equation is taken into account. At the macroscopic scale, this equation may be postulated to be the result of the inner product of Darcy's law and the seepage velocity. However, near the porous medium boundaries, this postulate seems questionable due to the spatial variations of the effective properties (velocity, permeability, porosity, etc.). In this work we derive the macroscopic mechanical energy equation using the method of volume averaging for the simple case of incompressible single-phase flow in porous media. Our analysis shows that the result of averaging the pore-scale version of the mechanical energy equation at the Darcy scale is not, in general, the expected product of Darcy's law and the seepage velocity. As a matter of fact, this result is only applicable in the bulk region of the porous medium and, in the derivation of this result, the properties of the permeability tensor are determinant. Furthermore, near the porous medium boundaries, a more novel version of the mechanical energy equation is obtained, which incorporates additional terms that take into account the rapid variations of structural properties taking place in this particular portion of the system. This analysis can be applied to multiphase and compressible flows in porous media and in many other multiscale systems.

I. INTRODUCTION

The mechanical energy equation is, together with the total mass and momentum equations, a key relationship for the study of fluid dynamics in many systems. As a matter of fact, the study of many engineering systems requires the use of the macroscopic versions of these three equations [ see Chapter 7 in Ref. [I]]. In the particular case of the macroscopic mechanical energy equation, a challenge usually lies on the prediction of the rate at which mechanical energy decreases in a system due to viscous dissipation. Interestingly, in natural porous systems, such as geological media, the information provided by the mechanical energy equation is usually not crucial, as, for instance, in the macroscopic study of groundwater flow. Therefore, the focus is directed to the macroscopic conservation equations of mass and momentum, i.e., the continuity equation and Darcy's law. This is due to the interest on the prediction of the piezometric head, which can be experimentally measured [see Chapter 4 in [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF]]. However, as emphasized by Celia and Norbotten [START_REF] Nordbotten | Geological Storage of CO 2 : Modeling Approaches for Large-Scale Simulation[END_REF] in engineering systems such as geothermal engineering, enhanced oil recovery, or geological CO 2 sequestration, the use of the continuity equation and Darcy's law is insufficient and it is necessary to take into account the information provided by the macroscopic mechanical energy equation, which is related to the total head [see Eq. 4.1.5 in Ref. [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF]]. In addition, in applications such as the disposal of high-level radioactive wastes, understanding of the thermohydromechanical processes taking place in the porous medium is crucial as remarked by Chen et al. [START_REF] Chen | [END_REF].

Furthermore, the study of the viscous dissipation term in the mechanical energy equation is also relevant in the study of nanofibrous media, in specific when the goal is to achieve outstanding properties such as stiffness, strength, and impact resistance [5].

The derivation of the macroscopic version of the mechani cal energy equation has been scarcely addressed in the litera ture. Some relevant works are those by Cushman [START_REF] Cushman | The Physics of Fluids in Hierarchical Porous Media: Angstroms to Miles[END_REF], who used the formalism of statistical mechanics to derive the average version of the total energy equation. A more recent derivation is available from the thermodynamically constrained averaging theory (TCAT) [START_REF] Gray | introduction to the Thermodynamically Constrained Averaging Theory for Porous Medium Systems[END_REF] by Gray and Miller, where this equation is used together with those from momentum and mass transport in order to derive a constrained entropy inequality, which is a key part of the TCAT theory. Following also a thermodynamic approach, Borja [8] analyzed the mechanical energy equation for unsaturated porous media. Furthermore, Zhu et al. [9], using the volume averaging method [IO], showed that the macroscopic form of the mechanical energy equation is equivalent to the result of taking the inner product of the seepage velocity with Darcy's law for single-phase unsteady flow in homogeneous porous media. In their analysis, they found that the term related to the rate of viscous dissipation is equal to Darcy's term dotted with the seepage velocity. Nev ertheless, this result may not necessarily hold near the porous medium boundaries, where rapid variations of the effective properties (velocity, permeability, porosity) take place as it is the case in many engineering applications. In addition, if more complicated situations are considered, such as multiphase flow or compressible flow, there are more reasons to doubt that the macroscopic energy equation in porous media is simply the result of the inner product of Darcy's law with the seepage velocity. Despite the plausible interest in these and other more complicated situations, the focus in this work is on the derivation of the mechanical energy equation using the method of volume averaging in the case of single-phase and steady state flow of an incompressible fluid that saturates a rigid porous medium. As a first approach, isothermal conditions are assumed, so only the mechanical part of the total energy equation is justified. Our aim is to first derive an expression for momentum transport that is sufficiently general so it can be applied everywhere in the system, i.e., in the bulk and near the porous medium boundaries, through which momentum exchange can occur (e.g., between a porous medium and a fluid or between two porous media). This expression can be subsequently simplified by imposing pertinent scaling postulates so the Darcy-Brinkman-Forchheimer (DBF) or just the Darcy-Forchheimer (DF) equations can be deduced from this general model. A relevant question to be addressed is if the result of upscaling the pore-scale mechanical energy equation is compatible with the result of taking the inner product of the DBF equation with the seepage velocity. This will allow us to determine up to which point the result from Zhu et al. [9] is valid and to provide a more general form.

The paper is organized as follows: first, we present the governing equations and boundary conditions for mass, momentum, and mechanical energy at the pore scale. The main steps involved in the volume averaging method are then recalled. Since our interest lies in the correspondence between the macroscopic mechanical energy and momentum transport equations (both near the boundaries and in the porous medium bulk) at the Darcy scale, we first direct the attention to the averaging of the momentum transport equation and then apply the volume averaging method to the mechanical energy equation. Finally, the corresponding discussions and conclusions are provided, taking into account extensions to other flow situations.

II. PORE-SCALE MODEL

Let us consider the steady, incompressible, and Newtonian flow of a single fluid-phase f3 saturating a rigid and homoge neous porous medium, as shown in Fig. I. Flow is assumed to take place in a range of the Reynolds number, justifying that both inertial and viscous forces play a significant role and the no-slip boundary condition is assumed to be applicable at the solid-fluid interface. On the basis of these starting assumptions, the governing equations for total mass and momentum in the f3-phase are V-v/J=O, (la)

PVtJ •VV fJ =-Vp fJ +PtJg+µV 2 v fJ . (lb)
The fluid velocity is subject to the interfacial boundary condition:

(le)

To complete the statement of the boundary-value problem, it is necessary to provide the boundary conditions at entrances and exits of the macroscopic system (d fJ ,e), which can be written as

V fJ = V fJ ,e, at dfJ,e, (ld)
where v fJ,e is assumed to be a known function of position. The solution of this boundary value problem provides the pressure and velocity fields everywhere in the {3-phase from which one could analyze the different terms in the pore-scale mechanical energy equation that is formally obtained after a dot product of the momentum equation with the velocity, yielding

�V • (v�v tJ ) = -V • (p fJ v fJ ) + pg• v fJ ] 2 2 T + 2 µ,V v fJ -µVv fJ : Vv/J. (2)
Here, as in the remainder of the development, the nested convention is adopted for the double inner products, i.e., A: B = AijBj; (= B: A). In addition, fluid density and viscosity are assumed to be constants.

Although achievable, the pointwise detailed solution on PfJ and v fJ may not be of interest while a macroscopic description is more relevant. As a consequence, there is a necessity for a systematic filtering of nonredundant information arising from the pore-scale model, a procedure that is referred to as upscaling [START_REF] Wood | [END_REF]. In this way, it is desirable to derive upscaled models starting from the pore-scale equations by means of a convenient averaging procedure. There are many methodologies to carry out this task such as the method of moments, homogenization, or, more recently, the ther modynamically constrained averaging theory. In this work, the method of volume averaging is employed to derive the upscaled models that are valid both at the porous medium bulk and near its boundaries. Some salient features of this method is that it not only provides the means to derive the upscaled models but also allows us to clearly identify the associated scaling postulates involved in the derivation of the models and incorporates a closure scheme to predict the associated macroscopic coefficients appearing in the model [12]. To have a clear view of the above, the main steps of the method of volume averaging are described in the following section.

III. OVERVIEW OF THE VOLUME AVERAGING METHOD

In Fig. 2 we show a scheme that describes the method of volume averaging. This scheme should be read from up to down and from left to right. In the latter case, the colors illustrate the systematic loss of redundant information, which we refer to as upscaling. Throughout this work, we will follow Wood [I I] and make a distinction between the mathematical operation of averaging and upscaling, because the first one does not require any reduction in the number of degrees of freedom involved in the model. For this reason, Fig. 2 shows two routes that can be followed to obtain different average models. The first route is followed by the continuous lines and it describes the classical method of volume averaging as used in Ref. [START_REF] Whitaker | The Method of Volume Averaging[END_REF], which leads to a closed upscaled model. The dashed lines in Fig. 2 indicate a second route that leads to an average equation that does not involve upscaling and is known as the one-domain approach [13], which can be further simplified to also yield a closed upscaled model as in the first route. In this work, we will follow this last approach. The rest of this section is devoted to explain in more detail each step represented in Fig. 2.

The first step of the method is the statement of the pore-scale model, which relies on a set of starting assumptions (number 1 in Fig. 2); this step was already described with the associated hypotheses in the previous section. Next, an averaging domain, Y (of norm V), such as the one shown in Fig. 1, is defined in such a way that it contains portions of all the phases involved in the analysis (i.e., Y = ytJ + % ) (letter a in Fig. 2). In terms of the averaging domain, the superficial and intrinsic averaging operators are introduced for a piecewise function defined everywhere in the ,8-phase, 1/1 {3 , as follows:

(1/l{J)lx = � f 1/l{Jlr dV, yE'1p(x) (3a) (3b)
Throughout this work, the equations written in terms of the piecewise function 1/1 {3 are applicable to scalar, vectorial, or tensorial quantities. In the above equations, the vector x locates the position of the centroid of the averaging domain whereas the position vectors y and r locate points in the ,8-phase relative to x and to an arbitrary reference system, respectively, as shown in Fig. 3. The superficial and intrinsic averaging operators are related by the Dupuit-Forchheimer equation: [START_REF] Chen | [END_REF] where .s(x) = V fJ (x)/ Vis the volume fraction of the ,8-phase contained in the averaging domain.

The superficial averaging operator is subsequently applied to the governing equations at the pore scale. This operation requires interchanging differentiation and integration (letter b in Fig. 2), which can be achieved by means of the general transport theorem [START_REF] Bird | Transport Phenomena[END_REF],

f n • 1/J fJ wdA, (Sa)
w being the velocity of S'i f3 a , and the spatial averaging theorem [14]:

(Vi/lfJ)lx = V(i/lfJ)lx + � f nijl fJ dA. (Sb) yEdp0(x)
For the case of the average of the divergence of a vectorial function, the spatial averaging theorem is slightly modified as indicated in Eq. (1.2-23) in Ref. [START_REF] Whitaker | The Method of Volume Averaging[END_REF]. It is worth mentioning that the interfacial integrals involved in the above theorems often allow direct substitution of the corresponding boundary conditions. At this point in the analysis, the average model involves integrals of pore-scale quantities and their averaged counterparts. We will refer to the integrals of pore-scale quantities as filters of information as suggested by Whitaker [START_REF] Whitaker | The Method of Volume Averaging[END_REF]. To eliminate the remaining pore-scale functions, the spatial decomposition introduced by Gray [15] is used (letter c in Fig. 2): fJ -

1/r{Jlr = (1/r fJ ) Ir+ 1/r{Jlr- (6)
For the remaining of this work, we will make a distinction between average quantities evaluated at r or at x, but we will no longer specify that pore-scale quantities (i.e., ifr fJ and 1frfJ) are evaluated at r. In this way, the resulting average model is unclosed because it involves two unknowns (the deviations and average quantities) and only one equation. To overcome this issue, it is necessary to derive and formally solve a boundary-value problem for the deviations fields, which we refer to as a closure problem. This is achieved on the basis of the spatial decomposition given by Eq. ( 6) and we will illustrate this in the following section. If no assumptions are made, then we can qualify this closure problem to be exact because it contains all the information from the pore scale. As explained by Wood and Valdes-Parada [12], despite the considerable complexity of this problem, it is possible to formally solve it using integral equation formulations based on Green's functions. If this solution is substituted into the filters of the unclosed model (letter d in Fig. 2), then the resulting average (but not upscaled) equation is the so-called one-domain approach because it can be applied everywhere in the system, i.e., in the bulk of the porous medium and near its boundaries. This completes the nonclassical route of application of the volume averaging method as explained above. Certainly, one may apply scaling postulates in order to obtain simplified versions of the one-domain approach and to even obtain the resulting simplified upscaled model that is only applicable in the bulk of the porous medium.

Analyzing this modeling approach, it is not hard to deduce that the kernel of the upscaling process relies on the systematic reduction of the information contained in the closure problem solution. Therefore, if reasonable scaling postulates (usually in the form of length-scale constraints and assumptions) are adopted (number 2 in Fig. 2), then one may significantly reduce the amount of information carried by the deviations variables and yet still capture the essential (i.e., nonredundant) information from the pore-scale model. The formal solution of this simplified closure problem can be substituted into the filters of the unclosed average model (letter d in Fig. 2) to obtain its closed form. However, the resulting model is still subject to additional simplifications that are based on complementary scaling postulates (number 3 in Fig. 2) in order to obtain a simplified closed upscaled model. In the following section, we will follow the methodology outlined here, i.e., we first derive the one-domain approach models for mass and momentum transport in porous media. Afterwards, pertinent scaling postulates are introduced in order to derive the corresponding models that are valid only in the bulk of the porous medium. Then a similar procedure is applied to the mechanical energy equation.

IV. UPSCALING MASS AND MOMENTUM TRANSPORT

A. Averaging

Let us, for the moment, direct the attention to mass transport and apply the superficial averaging operator as defined by Eq. (3a) to Eq. (la). The resulting equation is

(V • V fJ) Ix = 0. (7)
As explained in the previous section, with the aim of interchanging spatial differentiation and integration, we use the spatial averaging theorem in order to obtain:

V • (v fJ )lx = 0. (8)
Here we have taken into account the nonslip boundary condition that was imposed at the solid-fluid interface. It should be noted that this result is written in terms of the seepage velocity, which is the common form of averaging for the velocity found in the literature. As will be shown below, for the derivation of the closure problem, it is necessary to express the continuity equation in terms of the intrinsic averaged velocity. This is easily achieved by means of the Dupuit-Forchheimer relation [Eq. ( 4)]:

V • (v fJ ) /J lx = -V ln .s • (v fJ ) /J lx - ( 9 )
Directing the attention to momentum transport, let us apply the superficial averaging operator to Eq. (lb), taking into account the continuity equation, to obtain: p(V • V fJ V fJ )lx = -(V P fJ )lx -.spg + µ,(V 2 v fJ )lx - (10) The result of using the spatial averaging theorem to the transport terms in the above equation is pV

• (V fJ V fJ )lx = -V (p fJ )lx + .spg + µ,v' 2 (v fJ )lx + V 1 I n • ( -I PfJ + µ,Vv fJ )dA (11)
or, in terms of intrinsic averages, pV • ( c(V fJ V fJ ) /J Ix ) = -.sV ( PfJ j/J Ix+ .spg + µ,.sv'

2 (v{Jj/J Ix + µ,V • ( V.s(v fJ llx ) + � f n-[-l(p fJ -(P fJ ) /J lx ) yEd'fa(X)
With the aim of eliminating the pore-scale pressure and velocity contained in the surface and volume integrals (i.e., the filters) from this expression, we use the spatial decomposition defined in Eq. ( 6) so the result can be written as pV • ( .s((v fJ ) /J lr(V fJ ) /J lr + (v fJ ) /J lrV fJ + V fJ (V fJ ) /J lr + V fJ V fJ ) /J lx )

-.sV (p fJ ) /J lx + .spg + µ,.sv'

2 (v fJ ) /J lx + µ,V • ( V.s(v fJ ) /J lx ) + V (13)
It should be noted that each term in the above equation is evaluated at x. In addition, in the surface integral of the above equation, we used the following definition:

B. Exact closure problem

Equation ( 13) is the unclosed average model. Our goal now is to derive expressions for the velocity and pressure deviations. With this in mind, it is necessary to derive and formally solve the boundary-value problem that governs the fields of Pf! and v fJ • Substitution of Eq. ( 6) into the pore-scale continuity equation [Eq. (la)] yields

V-v fJ = -V • (v fJ ) fJ lr � volumetric (I Sa)
Directing the attention to the momentum transport equation, let us divide Eq. ( 13) by E, evaluate the resulting expression at x =r and subtract it to Eq. (lb) in order to obtain:

= -Vh + µ,V 2 v fJ -µ,E-1 V • (VE(V fJ ) fJ lr ) 1 / n-[-lh+µ,vvtJJdA vtJ zEd µ0 (r) I n • [-l�(PtJ) fJ + µ,V(�(v fJ ) fJ )JdA .
v fJ zEd µ0 (r)

volumetric source (lSb)
In this last result, we introduced the dummy vector z, in order to make a distinction with the position vector r in the volumetric and interfacial integrals.

The corresponding boundary conditions for the deviations, which result from substituting the spatial decomposition as given by Eq. ( 6 The formal solution of this problem is feasible using inte gral equations formulations based on Green's functions, as explained in Ref. [12], and it can be expressed as follows: In the above equations, G v and g P represent the corresponding Green's functions for the velocity and the pressure, respec tively, and r fJ represents a linear combination of the volumetric source terms identified in the differential equations (lSa) and (I Sb).

/ G v •

C. One-domain approach

Now that we have the solution of the closure problem, the average model given by Eq. ( 13) is complete. For the sake of simplicity in presentation, we introduce the following definition: f fJ (X) = -pV • (E((v fJ llrV fJ + V fJ (V fJ ) fJ lr + V fJ V fJ ) fJ lx )

+ V I + µ,V (v fJ + �(v fJ ) fJ )JdA, ( 17
)
which, once inserted into Eq. ( 13), yields the final form: pV • ((V fJ ) fJ lr(V fJ ) fJ lr)lx

= -EV (P fJ ) fJ Ix + Ep g + J,l V 2 (v fJ ) Ix -µ,VE• V(E-1 (v fJ )lx ) + f fJ (x). ( 18
)
The structure of this equation resembles its pore-scale counterpart, as given by Eq. (lb); however, it is important to clearly identify the physical meaning of each term in the above equation. The term in the left-hand side, although not written in the traditional form, is the macroscopic inertial term due to the macroscopic convective acceleration; the first term in the right-hand side corresponds to the macroscopic normal surface force due to the pressure gradient; the next one is due to volume forces in the fluid phase (gravity). The term µ, V 2 ( v .e ) I x is usually known in the porous media literature as the Brinkman correction term and it partially captures the macroscopic viscous forces. The term µ,Vs • V (c 1 ( v .e ) I) has been identified as the second Brinkman correction by Ochoa-Tapia and Whitaker [16]. It complements the mapping of the viscous forces by taking into account also the spatial variations of s. It is worth noting that this term is not present in the Navier-Stokes equations and it is thus a consequence of the averaging process. This characteristic is also shared by the term f .e in which all the contributions involving deviations, as defined in Eqs. ( 6) and ( 14), are gathered. As it will be shown later, all the terms involving spatial variations of the porosity and macroscopic velocity vanish in the porous medium bulk; however, the same is not true for f .e because the deviations of the averages from their pore-scale counterparts are still present in the bulk. As a matter of fact, only in single-phase systems can this term be safely discarded.

It should be noted that, up to this point, no scaling postulates have been introduced for the derivation of the average model, being consistent with the one-domain approach route depicted in Fig. 2. On the one hand, we have gained generality in the sense that the average model is applicable everywhere in the system; on the other hand, the complexity of this model may prevent its use in practical applications. Indeed, it appears to be simpler to perform direct numerical simulations from which the f .e term can be computed than using the solutions given in Eqs. (16). Motivated by this, we will systematically adopt scaling postulates with the aim of reducing the complexity of the one-domain approach model in the rest of this section.

D. Scaling postulates 1 and simplified closure problem

From the scheme illustrated in Fig. 2, a first set of scaling postulates can be imposed to simplify the closure problem (number 2 in Fig. 2). Therefore, as a first point of simplification, let us constrain the size of the averaging volume, r 0 , to obey the inequality e « ro « L, (19) where e represents the largest of the characteristic lengths associated to the pore scale while L denotes the smallest of the characteristic lengths associated to the macroscale. In this way, e = max(e .e ,e a-), where e .e and e a-denote the characteristic lengths associated to each phase as sketched in Fig. 1. It is worth mentioning that, using a phase-indicator function, it is possible to define additional relevant quantities such as the variance and the autocorrelation function, which allow determining the characteristic length scales in a more precise manner as explained by Wood and Valdes-Parada [12]. The interested reader is referred to this work for further details. Moreover, it is crucial pointing out that if there is no separation of length scales, then it is not possible to define a representative averaging domain and no further simplifications can be applied to the derivations. As explained by Quintard and Whitaker [17], for disordered media it is reasonable to assume that the size of the averaging domain must be several orders of magnitude larger than the characteristic length associated to the pore scale. However, there may be situations, such as fractal media, in which this assumption can fail. The remainder of this work is thus limited to systems in which it is reasonable to assume that the separation of length scales given in ( 19) is applicable. Under these conditions, we may refer to the averaging domain as a representative elementary volume (REV) [START_REF] Bear | Modeling Groundwater Flow and Contaminant Transport[END_REF].

Analyzing Eq. ( 18), taking into account the definitions given in Eq. ( 17), we observe that the one-domain approach model contains average quantities evaluated at r and at x. Actually, one may use Ta ylor-series expansions in order to express average quantities evaluated at r only in terms of average quantities evaluated at x [START_REF] Whitaker | The Method of Volume Averaging[END_REF], that is, (i/J .e ) .B lr =(i/J .e ) .B lx + y . V (i/J .e ) .B lx + ½YY : VV (i/; .e ) .e I x+••• .

(

) 20
Performing an order-of-magnitude analysis to the second and third terms on the right-hand side of this result leads to y • V(i/J .e ) .e lx = o (� (i/J .e ) .e lx ),

(1/r p )P yy: VV(i/; .e ) .e lx = o(L rJL (i/J .e ) .e lx ), (21b)

V(i/r p )P (1/r p )P where L; is the characteristic length associated to the spatial variations of the i quantity (i = (i/J .e ) .B , V (i/J .e ) .B ). To avoid oversimplifications for the problem under consideration, we associate the characteristic lengths Lv(,tr p )P and L(t/i p )P to the width of the transition zone near the porous medium boundary, say, O w ry • As determined by Valdes-Parada et al. [18], the width of this zone is roughly 20r 0 , indicating that there is at least one order of magnitude of difference between 8 w ry and r 0 . On this basis, we may accept the approximation (22) This implies that the average term on the right-hand side of the spatial decomposition given in Eq. ( 6) can now be evaluated at x. This allows us to consider average quantities as constants within (surface or volumetric) integrals. A direct consequence of this approximation is the following average constraint for the deviations fields:

Under these conditions, the closure problem can be signifi cantly reduced. The details of this simplification are given in Appendix, and it suffices here to present the simplified closure problem as

V • v .e = 0, (24a)
PV .e • Vv .e = -V f5.e + µ,V 2 v .e

I (24c)

At this point, we reduce the domain of solution of the closure problem to a periodic unit cell instead of the entire macroscopic domain. Actually, this simplification is more a convenience than a necessity and by no means constrains the use of the resulting average model to periodic systems, which are rarely encountered in practice. In this way, at the entrances and exits of the unit cell we can impose the following periodic boundary condition:

(24d)

where I; denotes the lattice vectors of the unit cell. Finally, it is worth recalling that the pressure deviations fields are bounded by the integral constraint (24e) Indeed, the velocity deviations satisfy a similar constraint, but such a constraint is unnecessary to have a well-posed problem in which the only remaining source term can be taken to be constant within the unit cell domain.

The formal solution of this problem can once again be obtained using integral equations formulations based on Green's functions, and for the sake of simplicity in notation, we express the solution as

V fJ = B fJ • (v fJ ) fJ l x , Pf3 = f,l b fJ • (v fJ ) fJ Ix- (25a) (25b)
Notice that this solution is a simplified version of the formal solution given in Eqs. ( 16); however, since now there are no volumetric sources and the boundary conditions at the entrances and exits are homogeneous, the only remaining source is the one located at the interface. In this way, the so-called closure variables B fJ and b fJ can be conceived as integrals of the corresponding Green's functions associated to the velocity and pressure, respectively. These closure variables solve the following boundary-value problem:

V • B fJ =0 !!... v fJ • V B fJ = -Vb fJ + V 2 B fJ , f,l B fJ = -I, at d fJ a, lf fJ (r) = lf fJ (r + I;); lf fJ = B fJ ,b fJ , ( b fJ ) fJ l x = 0. (26a) (26c) (26d) (26e)
Note that these equations correspond to Eqs. (2.22) in Ref. [19].

It is now opportune to return to the average momentum equation by first noticing that, under the assumptions adopted so far, f fJ can now be reduced to

f fJ (x) = -pV • (. s ( V fJ V fJ ) fJ l x ) + � f n-[-lh+µ,Vv
(v fJ ) fJ lx • J fJ • (v fJ ) fJ l x ) 8 = -V (p fJ ) fJ lx +pg+ t: v' 2 (V fJ ) Ix 8 -µV ln .s • V(.s-1 (v fJ )l x )-H i 1 • (v fJ )lx - (28)
In this result, we introduced the fourth-order tensor J fJ and the second-order tensor H fJ , which are defined in terms of the closure variables as follows:

J fJ = II+ (Br B fJ ) fJ l x ' (29a) -1 1 .sH = -- fJ v fJ f n • (-lb fJ + V B fJ )dA. (29b) yEJdp0 (x)
From a physical viewpoint, the fourth-order tensor J fJ helps grouping all the macroscopic inertial contributions in a single term. Furthermore, the term ( sr B fJ ) fJ I x is an inertial filter (in the sense of Whitaker [START_REF] Whitaker | The Method of Volume Averaging[END_REF]) that corrects the traditional inertial term by taking into account the contributions from the pore scale. In addition, the second-order tensor H fJ , which is a stress filter, may be regarded as an apparent permeability tensor that corresponds to a position-dependent Darcy-Forchheimer tensor. The introduction of these effective-medium coefficients is quite convenient because it allows us to separate the macroscopic inertial and viscous contributions in Eq. ( 28).

Before moving on, it should be noted that the closure problem given by Eqs. ( 26) is quite complicated because it involves an integrodifferential equation. With the aim of deriving a purely differential closure problem, we make the following changes of variables as suggested by Whitaker [19]:

E fJ = .s-1(s fJ +I)• H fJ , e fJ = .s-1 b fJ • H fJ . (30a) (30b)
These variables solve the following differential closure problem: From the above, it may appear that the closure problem requires the solution of the pore-scale problem in order to obtain the pointwise velocity that is present in the inertial term of Eq. (31b), making this boundary-value problem unclosed. However, it has been recently shown [20] that this is not the case, because one may use Eqs. ( 6) and (25a) in order to obtain a closed form of this problem.

V • E fJ = 0, !!... v fJ • VE fJ = -Ve fJ + V 2 E fJ + I,
In addition, from the constraint that the intrinsic average of the velocity deviations is zero, it follows that (3lf) This final form of the closure problem was recently solved by Lasseux et al. [21] for several flow angles and Reynolds num ber values on some model structures. Furthermore, following the derivations shown in Appendix from the study of gas-slip flow in porous media by Lasseux et al. [22], it is not hard to deduce that the H ,g tensor can be split into a symmetric and a skew-symmetric part, when restricted to the porous medium bulk (i.e., at x = x w ), as follows:

E. Scaling postulates 2 and closed upscaled model

With the simplifications to the closure problem adopted above, we have reached the reduced solution given by Eqs. (25) and, consequently, the effective-medium coefficients J ,g and H ,g can now be readily computed from the knowledge of closure variables fields. Despite the progress gained so far, the one-domain approach, as given in Eq. ( 28), still involves position-dependent macroscopic coefficients. Nevertheless, the scaling postulates introduced in the previous section have led to a considerable loss of information from the pore scale. Therefore, from the scheme in Fig. 2, we shall refer to this simplified version of the one-domain approach as a closed upscaled model with position-dependent coefficients. Actually, this simplified version of the one-domain approach has been used in the past for the derivation of jump conditions between a porous medium and a fluid under noninertial conditions [23]. However, in practice, it is more common to use effective-medium equations involving constant coef ficients and to account for this simplification by introducing appropriate jump conditions. This motivates the introduction of additional scaling postulates that are directed to further simplify the closed upscaled model and, as indicated in Fig. 2, we will refer to the resulting model as a simplified upscaled model.

With the above goal in mind, we first note that Eq. ( 28) involves three effective-medium coefficients, namely the volumetric fraction of the fluid-phase within the averaging volume, c(x), the apparent permeability tensor H ,g (x), and, finally, the tensor J ,g (x) that modifies the dyadic product of the average velocities in the inertial term. Outside the transition region that is located near the porous medium boundaries, i.e., in the porous medium bulk, these coefficients are constants and we denote the last position where this assumption is valid as Xw. In the transition region, we thus propose the following Taylor-series expansion of the coefficients about x w :

rp(x) = rp(x w ) + (x -Xw) • V rp(x w ) + ½ (x -Xw)(X -Xw) : VV <p(Xw) + • • • , ( 33
)
where <p = £, H ,g , J ,g . From this expansion, we note that if we truncate it at the first term, all the effective-medium coefficients in Eq. ( 28) can be replaced by their corresponding values at the porous medium bulk, which are position independent. Consequently, this equation reduces to the DBF equation with the inclusion of the inertial term:

pV • ( (v ,g ) ,'l l x • J ,g (X w ) • (v ,g Ji'l l x ) -V (p ,g ) ,'l Ix+ pg+ µV 2 (v ,g ji'l Ix -µ H fi 1 (x w ) • (v,g) l x (34)
Despite the drastic simplifications that we have imposed in order to obtain this equation, it may still be used near the porous medium boundaries as long as the errors induced by these approximations are compensated. One way to achieve this compensation is by taking more terms in the expansion given in Eq. ( 33) and another way is by the introduction of appropriate jump boundary conditions, which is the most common approach.

Finally, in the porous medium bulk, the velocity is position invariant and Eq. ( 34) is reduced to the form of the Darcy Forchheimer equation deduced by Whitaker [see Eqs. (v ,g )k = -H ,g (x w )

• ( V(p ,g ) ,'l l xwpg ). ( 35)

f,l
The main difference between the simplified upscaled models given in Eqs. (34) and ( 35) is that the first one may be used near the porous medium boundaries while the latter is constrained to the porous medium bulk. Before moving on, it is pertinent to point out that, even though the macroscopic viscous and inertial terms have been dropped, the pore-scale inertial and viscous forces contributions are still present in the above equation through the Darcy-Forchheimer term H ,g as it can be deduced from the analysis of Eqs. (31 ). Now that we have completed our analysis of the macro scopic momentum transport equation, we are in position to derive the average and upscaled versions of the mechanical energy equation, which is one of the main results of this article. This is the objective of the following section.

V. UPSCALING THE MECHANICAL ENERGY EQUATION

As mentioned in the Introduction, one of the objectives of this work is to derive a form of the macroscopic mechanical energy equation that is applicable both in the porous medium bulk and near its boundaries. This analysis raises interesting questions such as is the resulting expression equal to the inner product of the one-domain approach for momentum transport and the Darcy velocity? Is the result reported by Zhu et al. [9] in the bulk applicable? To this end, we will follow a similar structure of developments as the one used in the previous sec tion with the difference that we will not derive the one-domain approach version of the mechanical energy equation. Instead, we will first derive the closed upscaled version of this equation that is applicable at the porous medium boundaries and contrast it with the result of taking the inner product of the seepage velocity with Eq. (34). Finally, we will verify that when the simplified upscaled model for the mechanical energy equation because the upscaled model is simply the result of taking the inner product between the Darcy equation with the seepage velocity as shown in Eq. (50). This result may have not been achieved without knowledge of the symmetric and skew-symmetric parts of the apparent permeability tensor H tJ , as indicated in Eq. (32). Our analysis complements and verifies the result from Zhu et al. [ 9 ] .

As mentioned in the Introduction, the derivation in this work can surely be extended to other more complicated situations such as multiphase flow, compressible flow, and even to nonisothermal flow in porous media. As a final note, it is worth mentioning that the framework used in the present analysis, consisting of the derivation of a one-domain approach, which may be further simplified by the imposition of scaling postulates, is indeed extensible to other practical applications beyond porous media systems. Actually, the main message from this work is that the analysis in terms of energy balance near the boundaries of multiscale systems is incomplete if it deals only with the study of mass and momentum transport, because it is necessary to also consider the upscaled mechanical energy equation.

FIG. 1 .

 1 FIG. 1. Sketch of a porous medium including the averaging domain and the characteristic lengths.

  FIG.2. Diagram of the volume averaging method. The continu ous lines and arrows lead to the classical approach whereas the dashed lines and arrows lead to the one-domain approach.

  volumetr ic source + pV • [(v fJ ) fJ lrV fJ + V fJ (vtJllr + V fJ V fJ ] -E-1 P (V • [((v tJllz(vtJllz)lrJ volumetric source

  ) into Eqs. (le) and (ld), are V fJ = -(V fJ ) fJ lr , at d fJ rr, (lSc) V fJ = V fJ ,e -(vtJllr, at d fJ ,e • (lSd)

  r fJ dV y E'.f;,(x) Influence of volumetric sources / n • µ,V G v • (v fJ ) fJ l y dA y Ed µ 0 (X) Influence of the interfacial source + f n-µ,VG v•(V fJ ,e -(V fJ ) fJ l y )dA , (16a) y Ed p _,,(x) Influence of the entrances and exits Pf! = - f g P • r fJ dV y E'.f;,(x) Influence of volumetric sources / n • Vg p • (v fJ ) fJ l y dA y Ed p0 (x) Influence of the interfacial source + f n • Vg /J • (v fJ ,e -(v fJ ) fJ l y )dA . y Ed p _,.(x) Influence of the entrances and exits(16b)

  fJ ]dA. (27) yE.0fp0 (X) Substitution of the formal closure problem solution as given by Eqs. (25) into the above expression and subsequently into the momentum transport equation [Eq. (18)] yields !!...v • (.s

E

  fJ = 0, at tzl fJ a, lf fJ (r) = lf fJ (r + I;); lf fJ = E fJ ,e fJ , ( e fJ ) fJ Ix = 0.

  skew-symmetric contribution from inertial effectsIn the first term on the right-hand side of the above equation the superscript T3 permutes the first and third indices of a third-order tensor, i.e., A�i = Akji.

  (2.50) and (4.8) in Ref. 19]:

ACKNOWLEDGMENTS

C.T.P.G. is thankful to Universidad Aut6noma Metropoli tana for the scholarship received. F.J.V.P. is thankful to CNRS for the facilities provided to carry out a research internship with D.L. D.L. acknowledges financial support from CONACyT (Grant No. 25623) for a stay in UAM-Iztapalapa to work with F.J.V.P. during autumn 2016.

is used in the bulk, the result is the one reported by Zhu et al. [9].

Let us commence our derivations by applying the superficial averaging operator to Eq. ( 2) to obtain � (V • (v�v ,g )) l x = -(V • ( p ,g v ,g ) )lx +pg• (v ,g )lx + �µ(V 2 v�) l x -µ(vvr: Vv,,;) l x• ( 36)

Using the spatial averaging theorem and taking into account the nonslip boundary condition, we can interchange spatial integration and differentiation in order to express the above result as follows : To make further progress, let us adopt the scaling postulates I, so the formal closure problem solution given by Eqs. ( 25) is applicable. Therefore, taking into account the spatial decomposition in Eq . ( 6) and the changes of variables indicated in Eqs . ( 30), we can express the pore-scale velocity and pressure as follows :

v ,g = 0 ,g • (v ,g ) ,g

Ix,

where, for the sake of simplicity, we have introduced the following closure variables :

Using the definitions given in Eqs. ( 39), it can be easily deduced that

(vv;: Vv ,g ) ,g

l x = (v ,g ) ,g lx(V ,g ) ,g lx: ( ( VO ,,;

In the last equation, we have taken into account the separation of length scales that allows us to assume that V (v ,g ) iB where the physical meaning of each term has been clearly identified. As it was done with momentum transport, let us now adopt the second set of scaling postulates, that are focused to the upscaled model . From the previous section, we recall that the truncation of the expansion defined in Eq. ( 33), at the first term, allows us to approximate all the terms involving integrals of closure variables to their corresponding values in the bulk, i.e., at x w . Under these conditions, we may consider the volume fraction as a constant and divide both sides of the above equation by £ in order to obtain = -V. ( (p,g) ,g lx(v ,g ) ,g lx) +pg . (v ,g ) ,g lx -µV . ( (v ,g ) ,g lx. (d ,g O ,g ) ,g k . (v ,g ) ,g lx) + !µV 2 ((v ,g ) ,g lx(v ,g llx: (or. O ,g ) ,5 1 ) -µ(v ,g ) ,g lx(V ,g ) ,g lx: (( VO ,,; l

PAEZ-GARCIA, VALDES-PARADA, AND LASSEUX This equation is clearly not equal to the result of taking the inner product of Eq. ( 34) with the seepage velocity. Certainly, the above equation is the analog, in terms of mechanical energy, to the DBF equation, i.e., this equation can be used near the porous medium boundaries as along as there is an appropriate boundary condition that compensates the errors induced by the approximations involved in its derivation.

As a final point of analysis, it is necessary to verify if the upscaled model for the mechanical energy equation, when constrained to the porous medium bulk, does agree with Zhu et al. 's [9] result stating that the macroscopic mechanical energy balance is obtained from the dot product of the macroscopic momentum equation with the seepage velocity. With this in mind, let us fix the centroid of the averaging volume to values that lie in the porous medium bulk, i.e., x = x w , so the velocity can be assumed to be a constant, thus: 0 = -V(p ,'l ) ,'l lxw • (v ,'l ) ,'l lXm +pg• (v,'l) ,'l lxw -µ,(v ,'l ) ,'l lx,,, (v,'l) ,'l lx, , : ( ( VD ,'l l 3 : VD ,'l llx w , (

or, after a little rearrangement,

Directing the attention to the last integral term of the above equation and substituting the definitions given in Eqs. (4 0 ), one obtains:

or, using Eq. (32), ( ( VD ,'l f 3 : VD ,'l ) ,'l lx w cHtc(H�

(47) Substitution ofEq. (47) into Eq. ( 45) leads to 0 = (-v (P tJ lk, + pg -f,lc(v,'l) ,'l k •Ht)• (v,'l) ,'l lxw

where, for the sake of brevity in notation, we introduced:

As indicated in Eq. ( 32), the second-rank tensor (Er• [ V • ( vtJEtJ )] ) I is skew-symmetric and, since the dyad "' " s ,'l s ,'l is a symmetric second-order tensor, we can readily conclude that the last term on the right-hand side of Eq. ( 48) is zero. After multiplication by the porous medium porosity, this equation is simplified to which clearly is the inner product of Darcy's law with the seepage velocity, thus verifying the result from Zhu et al. [9] in the porous medium bulk.

VI. DISCUSSIONS AND CONCLUSIONS

In this work, we addressed the question about the correspondence between the upscaled version of the mechan ical energy equation and the inner product of the macroscopic momentum balance equation with the seepage velocity. The following points of discussion and conclusions are hence in order:

(i) Using the method of volume averaging, without the imposition of any length-scale constraints or assumptions, we derived an average model for momentum transport [Eq. ( 18)]. This one-domain approach, which is valid everywhere in the system (i.e., near the porous medium bulk and its boundaries), was the cornerstone for the derivation of upscaled models for momentum transport. By imposing a first set of scaling postulates, the one-domain approach gave rise to an upscaled model in which the effective-medium coefficients can be readily computed from the solution of ancillary closure problems in representative unit cells. Then a second set of more severe scaling postulates was imposed that lead to a version of the DBF equation that includes inertial contributions [Eq. ( 34)]. Finally, in the porous medium bulk, this last model reduces to the Darcy-Forchheimer model [Eq. ( 35)]. The approach for the derivation of these last two upscaled models differs from the classical applications of the volume averaging method reported by Whitaker [19].

(ii) Traditionally, Brinkman's correction to Darcy's law is regarded as an ad hoc addition to Darcy's law that allows the velocity to change with position near the porous medium boundaries, while keeping the effective medium coefficients constant. Actually, this is a contradiction, because, in general, the velocity and the effective-medium coefficients are both functions of position near the porous medium boundaries. This controversy can be solved by the inclusion of the second set of scaling postulates that give rise to the DBF model. Notice that, contrary to Whitaker [19], we have not imposed additional constraints dealing with the Reynolds number, so Eq. (34) maintains the macroscopic inertial contribution term.

(iii) The upscaled model for the mechanical energy equa tion resulting from applying the volume averaging method [Eq. ( 43)] was found not to be the result of taking the inner product of the DBF equation with the seepage velocity. However this expression of the macroscopic mechanical energy equation shares the same limitations (in terms of upscaling) as the DBF transport model. To the best of our knowledge, this complete form of the mechanical energy equation has not been previously reported in the literature. If this equation is to be used to compute the power required to reach a certain fluid flow near porous media boundaries through which momentum transport can occur, then it is necessary to acknowledge that the errors that are involved in its derivation should be compensated. There are at least two routes for achieving this goal: one way is to relax the scaling postulates involved in its derivation. For example, one may take more terms in the expansion of effective-medium properties defined in Eq. ( 33). An alternative is to derive the corresponding jump boundary conditions that incorporate coefficients that account for the information that was lost in the upscaled transport model.

(iv) Finally, in the porous medium bulk, it is not necessary to average the mechanical energy equation at the pore scale,

APPENDIX: SIMPLIFICATIONS OF THE CLOSURE PROBLEM

Our goal in this section is to detail the simplifications that are applicable to the closure problem associated to the one domain approach. On the basis of the separation of scales given in (19) and the corollaries presented in Sec. IV, the following simplifications are applicable to Eq. ( 15b):

(a) The last term on the left-hand side is reduced to

The last term on the right-hand side of this equation vanishes.

In this way, Eq. (15b) reduces to

To further simplify this equation, we use orders-of-magnitude estimates of the several terms that it contains. In performing this analysis, we take into account the interfacial boundary condition given by Eq. (15c), which indicates that the order of magnitude of the velocity deviations can be taken to be the same as the one for the intrinsic average velocity. Under these conditions, we have the following estimates:

where we have taken r 0 to be a reasonable estimate of the characteristic length scale of the spatial variations of the porosity. Since we have accepted the separation of length scales ,f, « r 0 , the above expression can be reduced to (A3a)

In a similar way, we have the following order-of-magnitude estimates:

In this way, on the basis of the length-scale constraints, l « ro; l « O w� , (A4) it is reasonable to assume that

Under these conditions, we can reduce Eq. (Al) to pV

Directing the attention to the left-hand side term of the above equation, we have From these estimates, we may assume that, due to the length scale constraint £ « r 0 , the interfacial source is much larger than the volumetric source. Under these conditions, the closure problem can be expressed as stated in Eqs. (24) in the main text.