Catherine O'sullivan
email: cath.osullivan@imperial.ac.uk

Chloé Arson

Benoît Coasne

A PERSPECTIVE ON DARCY'S LAW ACROSS THE SCALES: FROM PHYSICAL FOUNDATIONS TO PARTICULATE MECHANICS

This paper puts forward a perspective or opinion that we can demonstrate Darcy's law is valid at any scale where fluid can be modelled/analyzed as a continuum. Darcy's law describes the flow of a fluid through a porous medium by a linear relationship between the flow rate and the pore pressure gradient through the permeability tensor. We show that such a linear relationship can be established at any scale, so long as the permeability tensor is expressed as a function of adequate parameters that describe the pore space geometry, fluid properties and physical phenomena. Analytical models at pore scale provide essential information on the key variables that permeability depends on under different flow regimes. Upscaling techniques based on the Lippman-Schwinger equation, pore network models or Eshelby's homogenization theory make it possible to predict fluid flow beyond the pore scale. One of the key challenges to validate these techniques is to characterize microstructure and measure transport properties at multiple scales. Recent developments in imaging, multi-scale modeling and advanced computing offer new possibilities to address some of these challenges.

• 𝐴 cross sectional area • MD molecular dynamics

• 𝑛 porosity

• NMR nuclear magnetic resonance

• 𝑁 𝑝 number of adjacent pores

• 𝑃 pressure

• P average pressure

• PFGNMR pulsed field gradient nuclear magnetic resonance

• Q flow rate
• Q 𝑎𝑏 flow rate from pore 𝑎 to pore 𝑏

• QENS quasielastic neutron scattering

• 𝑟 radial direction origin at centre of cylinder

• 𝑅𝑒 Reynolds number

• REV Representative Elementary Volume

• 𝑠 specific surface area of system

• SPH smooth particle hydrodynamics inducing the flow [START_REF] Darcy | Les fontaines publiques de la ville de Dijon : exposition et application des principes à suivre et des formules à employer dans les questions de distribution d'eau[END_REF]. When Darcy's law is used in civil engineering applications, the flow rate is usually expressed in terms of a gradient of the total head, i.e. as:

v = - k 𝜂 • (∇𝑃 -𝜌g) (1)
where k is the permeability (sometimes termed the intrinsic permeability), ∇𝑃 is the gradient in pressure, the vector g represents acceleration due to gravity, 𝜌 is the fluid density, and 𝜂 is the fluid viscosity under given pressure/temperature conditions. In Eq. 1, the parameter k relates two vector quantities and so it is a tensor denoted 𝑘 𝑖 𝑗 . The unit of (intrinsic) permeability is the darcy, 1 darcy= 9.869233 × 10 -13 𝑚 2 or 0.9869233𝜇𝑚 2 . In geotechnical engineering, where Darcy's law is expressed in terms of a gradient in the total head (energy/unit weight) and the permeating fluid is usually water the term permeability (with units of 𝑚/𝑠) is generally used for K ′ = k𝜌𝑔/𝜂 [START_REF] Cashman | Groundwater Lowering in Construction A Practical Guide to Dewatering[END_REF]; arguably for clarity and precision the term hydraulic conductivity should be used for K ′ when water is the permeating fluid. The term permeance is used to designate the quantity K = k/𝜂, which is not restricted to water as the permeating fluid. Typical values of K ′ for geomaterials are given in Table 1.

The simple, linear relationship that is Darcy's law, is applied in physics and chemistry, materials and Earth sciences, mechanical and chemical engineering, soil and agricultural sciences, etc. In civil engineering, Darcy's Law is routinely used in geotechnical engineering design and analysis, including but not limited to the design of dewatering systems to control ground water during construction [START_REF] Cashman | Groundwater Lowering in Construction A Practical Guide to Dewatering[END_REF] and seepage analysis to design and assess levee systems [START_REF] Smith | Chapter 9 design[END_REF]). Darcy's law is also applied when interpreting well pumping tests to assess hydraulic conductivity in situ. Even though Darcy's law was historically introduced empirically 150 years ago, it was fully justified within the theoretical framework of homogenization at the 4 O'Sullivan, June 1, 2022

scale of continua at the end of the twentieth century (e.g., [START_REF] Boutin | Homogenization of coupled phenomena in heterogenous media[END_REF]). However, with the emergence of nanoscience, many fundamental questions have been raised regarding the validity of Darcy's law at vanishing length scales. In particular, considering that the flowing fluid can no longer be treated as a continuum medium at the nanometer scale, many poorly known phenomena become important as the molecular nature of the fluid comes into play [START_REF] Bocquet | Nanofluidics, from bulk to interfaces[END_REF].

At the same time, questions about Darcy's law have been put forward in the civil engineering literature. For example, [START_REF] Hansbo | Consolidation equation valid for both darcian and non-darcian flow[END_REF] has argued that Darcy's law is not valid in clay at low hydraulic gradients. [START_REF] Ejezie | Permeation behaviour of phpa polymer fluids in sand[END_REF] have highlighted the challenges associated with applying Darcy's law to predict the behaviour of non-Newtonian polymer fluids in sand, an application that is important to develop effective support for deep excavations.

Recent developments in numerical analysis and imaging, alongside an emerging emphasis on multi-scale analysis, allow us to re-examine Darcy's law. In this paper, we explore how we can bridge In this contribution we explore the general applicability of Darcy's law by considering flow in a single, cylindrical pore for different scenarios, including non-laminar flow and non-advective flow. This abstract model provides us with a clear basis to understand some of the limitations of the applicability of Darcy's law and we take it a basis to reflect on how to extend Darcy's law to poromechanical problems that involve strong couplings with chemical and physical processes that occur at the pore scale and below. We distinguish between analysis of fluid at the molecular and O'Sullivan, June 1, 2022 continuum scales and show that, once water can be modelled as a continuum, Darcy's law is valid at all scales so long as the permeability is expressed as the correct function of geometric, chemical, physical and mechanical parameters that fluid flow depends on.

VALIDITY OF DARCY'S LAW AT MULTIPLE SCALES

Considerations of scale must be addressed in even the most basic analysis of Darcy's law. The term v in Eq. 1 refers to a flow rate, Q, divided by a cross sectional area, 𝐴, orthogonal to the direction of flow, i.e. 𝑣 𝑖 = 𝑄 𝑖 /𝐴. The term v is sometimes referred to as the Darcy flux or Darcy velocity [START_REF] Cashman | Groundwater Lowering in Construction A Practical Guide to Dewatering[END_REF] or the discharge velocity [START_REF] Cedergren | Seepage, Drainage, and Flow Nets[END_REF]. The area, 𝐴, is occupied by both pore space and solid space, and so, clearly, without considering any local or global heterogeneities, the term v in Darcy's law is not the velocity of fluid flow. [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF] and [START_REF] Cashman | Groundwater Lowering in Construction A Practical Guide to Dewatering[END_REF] refer to v as an apparent velocity; v can also be seen as "a statistical summary of viscous flow" [START_REF] Ineson | Darcy's law and the evaluation of permeability[END_REF]. The actual pore water flow velocity will always differ from v. The mean pore velocity is given by v 𝑠 = v/𝑛 (where 𝑛 is the porosity)

and is termed the average linear groundwater velocity [START_REF] Cashman | Groundwater Lowering in Construction A Practical Guide to Dewatering[END_REF] or the seepage velocity [START_REF] Cedergren | Seepage, Drainage and Flow Nets[END_REF].

Sub-pore scale -numerical models

Modern numerical and experimental techniques enable us to consider the details of fluid flow in the pore space. At the molecular level, the arrangement of water in the pore space can be highly complex due to mechanisms including multi-layer adsorption, surface hydration and cation hydration (e.g. [START_REF] Zhang | What is the range of soil water density? critical reviews with a unified model[END_REF]; [START_REF] Coasne | Multiscale adsorption and transport in hierarchical porous materials[END_REF]). As illustrated schematically in Figure 1 (a), water in the pore space can be studied by molecular models using either Molecular Dynamics or Monte Carlo simulations with an appropriate water model (e.g. [START_REF] Pathirannahalage | Systematic comparison of the structural and dynamic properties of commonly used water models for molecular dynamics simulations[END_REF]).

The insight provided by considering these molecular mechanisms can inform our understanding of fluid flow in porous media. The exact dimension/resolution at which the continuum assumption applies depends on the fluid in question, the temperature and the pressure; for any fluid the continuum-discrete limit is approximately between 1nm and a few nm. As discussed in (Bocquet 6 O'Sullivan, June 1, 2022

and Charlaix, 2010), this limit can be determined from the correlation time of the stress-strain correlation function.

When the fluid is viewed as a continuum, direct simulation of fluid flow in the pore space can be achieved if the 3D topology of the pore space is available. This 3D topology data can be obtained from numerical simulations using the discrete element method (e.g. [START_REF] Knight | Computing drag and interactions between fluid and polydisperse particles in saturated granular materials[END_REF] . Whether the LBM or the FVM method is used, assumptions must be made regarding the boundary condition that is applied at the fluid-particle interface; typically a no-slip boundary condition is applied. In this context, it is important to realize that these methods implicitly or explicitly assume some specific regimes which often fail to capture the complexity of fluid dynamics in the vicinity of the surface (slippage, adsorption, non-viscous effects including memory effects, etc.). Some of these aspects will be covered in the next section. An additional key challenge in performing continuum simulations such as FVM or LBM is meshing the often complex topology of the pore space. The computational cost of the simulations restricts considerations to relatively small sub-volumes, however experimental work has shown that despite the resultant boundary effects, this type of model can capture the details of fluid flow with reasonable accuracy [START_REF] Sanvitale | Particle-scale observation of seepage flow in granular soils using piv and cfd[END_REF]. Results from these simulations can inform our understanding of Darcy's law at the macro-scale. While Darcy's law does not emerge from these models, they allow us to estimate the permeabilities (i.e. the proportionality constant between the flow and the pressure drop inducing the flow) for different thermodynamic and hydrodynamic regimes. In particular, by post-processing the fluid velocity field, one can calculate the overall permeability of an assembly of particles with a complex pore topology (e.g. [START_REF] Taylor | Sub-particle-scale investigation of seepage in sands[END_REF]) or the conductance that governs 7 O'Sullivan, June 1, 2022 flow in an individual pore throat (e.g. [START_REF] Sufian | Ability of a pore network model to predict fluid flow and drag in saturated granular materials[END_REF]). As discussed below, homogenization approaches can be applied.

Sub-pore scale -analytical models

While useful insight can be gained from numerical models that operate at the sub-pore scale, analytical models have the advantage of developing closed-form expressions where the relationships between key variables are explicit. By assuming that the pressure gradient inducing transport is small enough for the flow to be laminar, viscous and incompressible, Poiseuille-Hagen flow equations can be applied in media with simple pore network shapes (morphology/topology).

In what is perhaps the simplest scenario, we consider a cylindrical pore of length 𝐿 and diameter 𝐷 filled with an incompressible fluid having a (constant) viscosity 𝜂 and subjected to a pressure gradient ∇𝑃. For small Reynolds numbers, the flow can be assumed to be laminar as inertial contributions v • ∇v (where v is the velocity field) are negligible. In this case, the Navier-Stokes equation for momentum conservation reads [START_REF] Bird | Transport Phenomena[END_REF]:

𝜌 𝜕v 𝜕𝑡 = -∇𝑃 + 𝜂∇ 2 v (2)
For steady-state (stationary) conditions where 𝜕v/𝜕𝑡 = 0, the equation reduces to ∇𝑃 = 𝜂∇ 2 v.

Taking 𝑥 to be the direction of the pressure gradient, i.e. ∇𝑃 𝑥 = -Δ𝑃/𝐿, and 𝑟 as the radial direction with origin at the centre of the cylinder, one obtains a velocity field v = 𝑣 𝑥 (𝑟)e x (where e x is a unit vector in the 𝑥 direction). Integrating twice with respect to 𝑥 leads to 𝑣 𝑥 (𝑟) = -Δ𝑃𝑟 2 /2𝐿𝜂 + 𝐶 where the constant 𝐶 is calculated from the velocity boundary condition imposed by solid surface.

Macroscopically, this boundary condition at the pore surface usually corresponds to a "no slip" condition so that 𝑣 𝑥 (𝑟 = 𝐷/2) = 0 and hence 𝐶 = Δ𝑃/8𝐿𝜂 × 𝐷 2 (as discussed in the previous section, this no-slip boundary condition is usually assumed to be valid, without further justification, in mesoscopic approaches such as FVM and LBM). In turn, this leads to the following velocity profile for a cylindrical pore:

𝑣 𝑥 (𝑟) = Δ𝑃𝐷 2 16𝐿𝜂 1 - 2𝑟 𝐷 2 (3) 8 O'Sullivan, June 1, 2022
Hence, as illustrated in Figure 2, at any given position, 𝑟, there is a linear dependency of the velocity 𝑣 𝑥 (𝑟) upon the pressure gradient, Δ𝑃/𝐿. The latter equation allows deriving the average velocity over the whole cross-section area (𝑣 = Δ𝑃𝐷 2 /32𝜂𝐿) and the maximum velocity in the pore center (𝑣 max = Δ𝑃𝐷 2 /16𝜂𝐿).

The velocity boundary condition imposed by the solid surface depends on the detailed fluid/solid wetting properties and the structure of the solid surface. Under significant wetting, the large fluid/solid interactions lead to no-slip conditions and the velocity vanishes at the pore surface, leading to Eq. 3. In contrast, for partial wetting, the fluid slips at the pore surface with a velocity field that does not vanish when 2𝑟 → 𝐷. In this case, the so-called slip length 𝑏 corresponds to the distance at which the linearly extrapolated velocity vanishes, i.e. -𝑏𝜕𝑣 𝑥 /𝜕𝑟 = 𝑣 𝑥 (see Figure 2.b over the length of 𝑏). When this boundary condition is applied, the constant obtained following double integration becomes 𝐶 = Δ𝑃𝐷 2 /16𝐿𝜂 × [1 + 4𝑏/𝐷], giving the following velocity profile:

𝑣 𝑥 (𝑟) = Δ𝑃𝐷 2 16𝐿𝜂 1 + 4𝑏 𝐷 - 2𝑟 𝐷 2 (4)
When compressible fluids are considered (e.g. gases, supercritical fluids, etc.), special attention should be paid as to whether or not the flow can be treated as incompressible. Flow can be classified as compressible or incompressible by considering the dimensionless Mach number, defined as Ma = 𝑣/𝑐, where 𝑐 is the sound velocity in the fluid and 𝑣 can be taken as the maximum velocity.

For Ma < 0.3, the fluid flow can be regarded as incompressible because the flow is not affected by the pressure wave induced upon fluid displacement. On the other hand, for Ma > 0.3, fluid compressibility affects the flow and the Hagen-Poiseuille flow provided above must be modified.

Compressible flow behaves with an apparent slip length, 𝑏, defined as 𝑣 𝑥 = -𝑏𝜕𝑣 𝑥 /𝜕𝑟, as before, however the slip length 𝑏 scales with the inverse of the pressure, i.e. 𝑏 = 𝑏 0 /𝑃, where 𝑏 0 is a system-dependent constant while 𝑃 is the average pressure between the upstream and downstream pressures. The velocity profile for such a compressible flow in a cylindrical pore of diameter 𝐷 can O'Sullivan, June 1, 2022 then be expressed as [START_REF] Bird | Transport Phenomena[END_REF]:

𝑣 𝑥 (𝑟) = Δ𝑃𝐷 2 16𝐿𝜂 1 + 4𝑏 0 𝐷𝑃 - 2𝑟 𝐷 2 (5)
This equation corresponds to the so-called Klinkenberg effect which can describe flow that is an intermediate regime between Knudsen diffusion and viscous flow. As expected, for large pressures 𝑃, the fluid becomes incompressible so that slippage effects captured in 𝑏 0 /𝑃 vanish and Eq. 5

becomes equivalent to the Poiseuille flow predicted under the no-slip hypothesis (Eq. 3).

Both Eq. 3 and 4 predict a linear relationship between 𝑣 𝑥 (𝑟) and Δ𝑃/𝐿 for incompressible flow whether or not there is slip at the fluid-solid interface. While it may appear that Eq. 5 again gives a linear relationship between 𝑣 𝑥 (𝑟) and Δ𝑃/𝐿 for incompressible flow, the situation is more complex as the constant of proportionality depends upon 𝑃 as well as upon the pressure variations 𝑃(𝑥)

along the pore length 𝐿. However, in most situations, provided the pressure variations are smooth, one can write 𝑃(𝑥) ≃ 𝑃 + 𝑥∇𝑃 where ∇𝑃 is essentially constant so that the proportionality factor in Eq. 5 only depends on 𝑃 and ∇𝑃.

Upscaling from the Hagen-Poiseuille equations

The Hagen-Poiseuille equations above were derived for a simple pore geometry so that a welldefined velocity boundary condition can be imposed. When considering more realistic materials, a more complex porous volume geometry must be used. Using the same linear response regime formalism and assuming the fluid flow to be incompressible and purely advective, the Hagen-Poiseuille equation can be extended to any real medium.

For ideal pore geometries such as slit and cylindrical pores, the permeability k can be inferred by mapping Darcy's law onto the corresponding Poiseuille flow. For a porous medium with a porosity 𝑛, the flux J = 𝜌𝒗 is given by:

J = - 𝜌𝑛k 𝜂𝜏 • ∇𝑃 (6)
where the dimensionless factor 𝜏 is coined as tortuosity. This parameter ensures that Eq. 6 holds 10 O'Sullivan, June 1, 2022 when relating the observed flow J to the other physical constants that appear in the equation.

In practice, this tortuosity correction accounts for several effects that pertain either to the pore network geometry (e.g. flow paths that depart from a straight line, non-cylindrical pore shapes) or adsorption and surface phenomena (e.g. complex pore surface hydrodynamic boundary, residence times).

A particular form of Darcy's law is obtained when using the Kozeny-Carman relation for 1D

flow in which permeability is a scalar noted 𝑘 which is a function of the geometry of the porous space. Using the specific surface area of the system 𝑠, the Kozeny-Carman relation writes:

𝑘 = 1 5𝑠 2 × 𝑛 3 1 -𝑛 2 ∼ 𝐷 2 180 × 𝑛 3 1 -𝑛 2 (7)
where the second equality corresponds to the situation of spherical pores having a size 𝐷 so that 𝑠 = 6/𝐷. The Kozeny-Carman expression of permeability is used when applying Darcy's law to reflect the effect of porosity, pore size and pore shape (tortuosity) on permeability. In other words, the Kozeny-Carman model gives an estimate of the effective permeability based on simple descriptors of the microstructure. Thus, it is not generalizable to all porous media.

WHEN IS FLUID FLOW NOT A LINEAR FUNCTION OF THE GRADIENT OF PRESSURE?

Darcy's law is a linear response equation which can be seen as a specific expression of Onsager's phenomenological theory of transport applied to fluid flow in a porous medium [START_REF] Onsager | Reciprocity relations in irreversible thermodynamics. i[END_REF][START_REF] Barrat | Basic Concepts for Simple and Complex Liquids[END_REF]). Darcy's linear equation has proven to accurately describe laminar flow at several scales. This section reviews cases when the relationship between fluid flow is not a linear function of the gradient of pressure, including well-known problems of turbulent flow, gas flow, and fluid slippage along the pore boundaries.

Non laminar flow

It is expected that Darcy's law will apply only for a small Reynolds number Re (Re < 10, where 1D flow for the sake of simplicity but the description can be extended to 3D anisotropic systems:

-∇𝑃 ∼ 𝑎 ′ 𝑣 + 𝑎 ′′ 𝑣 2 (8)
where 𝑎 ′ and 𝑎 ′′ are constants. Physically, Eq. 8 accounts for non-linear effects as turbulence arises upon increasing Re. When such effects are included in Darcy's law, one obtains the so-called Darcy-Forchheimer law [START_REF] Nield | Convection in Porous Media[END_REF]:

- 𝑘 𝜂 • ∇𝑃 = 𝑣(1 + 𝑏 ′ 𝑣) = 𝑣(1 + 𝑏 ′′ Re) (9)
where the second equality on the right-hand side is derived by using Re = 𝜌𝑣𝐷/𝜂 (𝑏 ′ and 𝑏 ′′ are constants).

Transition from advective to diffusive flow

The assumption of advective flow is central to the derivation of Darcy's law. Depending on the flowing fluid -i.e. liquid, gas, supercritical fluid, the flow can be incompressible or compressible with underlying transport mechanisms of varying complexity that strongly depend on the applied thermodynamic boundary limits. In many situations relevant to practical applications, the flow rate complexity is a consequence of the dependence of the local density at a given position in the porous medium on the local pressure.

Knudsen diffusion corresponds to situations where fluid transport mostly occurs through diffusion of molecules colliding with pore walls as the fluid mean free path far exceeds the pore size.

In general, fluid transport in porous media can be rationalized by considering its mean free path 𝜆 ∼ 1/𝜌𝜎 2 where 𝜌 and 𝜎 are the fluid density and kinetic diameter, respectively (typically, 𝜎 is a fluid-specific quantity related to the scattering cross-section, which represents the probability that a molecule collides with another). For compressible phases such as gases or supercritical fluids, 𝜆 depends strongly on the applied pressure with values that can be either close to the pore size 𝐷 ("tight gas") or very large compared to 𝐷 (infinite dilution). The Knudsen number Kn = 𝜆/𝐷 is 12 O'Sullivan, June 1, 2022 the appropriate dimensionless number to describe the flow type as a function of 𝜆. Depending on the Knudsen number characterizing the flow, different regimes that deviate from Darcy's law can be encountered, as summarized here [START_REF] Ziarani | Knudsen's Permeability Correction for Tight Porous Media[END_REF][START_REF] Bear | Modeling Phenomena of Flow and Transport in Porous Media[END_REF]:

Viscous flow, Kn< 0.01. In this asymptotic regime, the flow is advective as described using Darcy's law. In this case, the flowing fluid is in the laminar regime (i.e. low Re).

Slip flow, Kn∼ [0.01 -0.1]. In this regime, the flow obeys Darcy's law corrected for Klinkenberg effect. As discussed above, this correction corresponds to compressibility-induced slippage at the solid surface.

Transition flow, . This transition regime corresponds to complex transport mechanisms where advection and diffusion are combined. Under such regimes, these two contributions sum up, as advection and diffusion occur in parallel [START_REF] Kärger | Diffusion in Nanoporous Materials[END_REF]. Upon increasing Kn, the flow evolves in a continuous fashion from purely advective to purely diffusive.

Knudsen diffusion, Kn> 10. In this asymptotic domain, the fluid mean free path largely exceeds the pore size so that Knudsen diffusion prevails with an underlying transport mechanism governed by molecular collisions with the solid surface.

Deviation from Darcy's law at the nanoscale

Adsorption effects can drastically affect the permeability and lead to a complex dependency of k on fluid density i.e. k(𝜌). Where this arises, upon imposing a pressure gradient across a porous medium to induce fluid transport, non linear effects are expected. This departure from Darcy's law stems from the strong adsorption at the pore surface. As explained in [START_REF] Falk | Subcontinuum mass transport of condensed hydrocarbons in nanoporous media[END_REF], these molecular effects make the viscous flow approximation at the heart of Darcy's law invalid. Without going too much into details, for fluid flow in a porous medium, Stokes' equation for momentum conservation in a laminar flow writes: Assuming for the sake of simplicity isotropic flow, the general solution of this differential equation is given by:

𝜌 𝜕v 𝜕𝑡 = -∇𝑃 + 𝜂∇ 2 v -𝜉v (
⟨ 𝑗 𝑥 (𝑘, 𝑡) 𝑗 𝑥 (-𝑘, 0)⟩ ∼ exp -𝜂𝑘 2 𝜌 -𝜉 𝑡 (11)
We note that the latter equation does not involve the pressure gradient because such time correlation functions are probed in the direction normal to the flow (i.e. perpendicular to the pressure gradient).

As illustrated in [START_REF] Falk | Subcontinuum mass transport of condensed hydrocarbons in nanoporous media[END_REF], in some cases, velocity fluctuations for nanoconfined fluids do not follow the behavior expected according to Eq. 11. This indicates that in such extreme confinement, the flow is not simply viscous as assumed in Darcy's equation. In this case, memory effects as defined in generalized hydrodynamics have to be taken into account. as Green-Kubo formula for the viscosity, can be derived in the framework of the linear response theory [START_REF] Barrat | Basic Concepts for Simple and Complex Liquids[END_REF]. In the specific context of viscosity, this equation indicates that the viscosity characterizes how fast the stress relaxes in the fluid particle (here, "particle" is understood in the context of hydrodynamics so that it corresponds to a small fluid element/volume in which the fluid can be treated as a continuum medium). Typically, the relaxation time 𝜏 𝑅 for the stress tensor is of the order of 10 -12 s even if the value can vary broadly depending on fluid molecular nature, temperature, density, etc. On the other hand, momentum transfer occurs over a characteristic time given by the decay time in the correlation function [START_REF] Barrat | Basic Concepts for Simple and Complex Liquids[END_REF]:

⟨ 𝑗 𝑥 (𝑘, 𝑡) 𝑗 𝑥 (-𝑘, 0)⟩ ∼ exp[-𝑡/𝜏] with 𝜏 = (𝑘 2 𝜌/𝜂 -𝜉) -1 . In the continuum (macroscopic) limit, 14 O'Sullivan, June 1, 2022
there is a marked time scale separation between these different phenomena as 𝑘 → 0, 𝜏 𝑅 ≪ 𝜏. In contrast, in extreme confinement and/or with strong surface adsorption effects, momentum transfers occur with wave vectors 𝑘 ∼ 1/𝐷 where 𝐷 is the pore size. As explained in [START_REF] Bocquet | Nanofluidics, from bulk to interfaces[END_REF], for such ultra-confined fluids, time scale separation is no longer valid as 𝜏 𝑅 ∼ 𝜏.

From a practical viewpoint, the viscous breakdown does not render Darcy's law invalid in the sense that proportionality between the flow rate and gradient pressure is still observed. However, strong molecular effects such as adsorption lead to permeabilities that are no longer an intrinsic property of the host porous medium. In particular, the permeability becomes a complex quantity that depends on both the solid and fluid properties as well as on the details of their interaction.

Permeability and Permeance

While Darcy's law can almost always be used empirically, the challenge in expressing permeability as a function of solid and fluid properties raises questions regarding the concept of permeability as a robust physical descriptor for flow in complex porous media. There are essentially two scenarios where the assumptions of viscous flow and ideal pore geometry become questionable: (1) For disordered pore morphology and pore network topology, regardless of the pore size, the permeability can no longer be written as a simple function of the pore size as the complex geometry does not make it possible to write well-defined boundary conditions. (2) In extreme confinement, such as in nanoporous media, the complex thermodynamics of the fluid in the vicinity of the solid surface is strongly coupled with the hydrodynamic flow. The permeability becomes a complex function of the fluid local properties (e.g., density, adsorption, dynamics).

In very small pores, where the fluid molecular granularity can no longer be ignored, additional questions arise regarding, for instance, the validity of macroscopic concepts such as viscosity.

While it is often assumed that 𝑘 only depends on the porous solid properties (typically, the pore size 𝐷 and/or specific surface area 𝑠), there is significant experimental evidence for complex permeability changes upon varying the fluid, density, temperature, etc. These effects were discussed by [START_REF] Coasne | Multiscale adsorption and transport in hierarchical porous materials[END_REF]. Due to the intrinsic complexity of these combined effects, additional phenomena occur at different scales (e.g. adsorption-induced poromechanical effects, memory effects in fluid flow, thermodynamical and dynamical phase transitions).

As a result, attempts to use theories including those relying on computational fluid dynamics tools may fail to establish robust expressions for the fluid permeability in a given porous medium. Combined geometric, chemical, physical and mechanical effects started to receive full attention with the advent of nanoporous materials (especially in the context of nanofluidics), but they are expected to apply to almost all practical situations as they are intrinsically linked to the thermodynamics and dynamics of confined and/or interfacial fluids.

As mentioned above, Darcy's law can be expressed using a single parameter K ′ = k𝜌𝑔/𝜂 known as hydraulic conductivity or K = k/𝜂 known as permeance. The expression K ′ = k𝜌𝑔/𝜂 implies that permeability k is an intrinsic property of the host porous material and that the use of a single hydraulic conductivity parameter K ′ couples the intrinsic properties of the fluid and solid. The same is true of the permeance K = k/𝜂 defined from momentum balance equations.

For example, using molecular modeling, [START_REF] Falk | Subcontinuum mass transport of condensed hydrocarbons in nanoporous media[END_REF] and later [START_REF] Obliger | Free Volume Theory of Hydrocarbon Mixture Transport in Nanoporous Materials[END_REF] showed that the 1D flow of different hydrocarbons (density, chain length) in disordered nanoporous solids cannot be reconciled with Darcy's law expressed as 𝑣 𝑥 /∇𝑃 = 𝑘/𝜂. On the other hand, the flow of confined molecules was found to be linearly related to the pressure gradient i.e. 𝑣 ∼ -𝐾∇𝑃.

The permeance 𝐾 is a complex function of the fluid density and molecule type. Yet, as expected for an incompressible fluid, the permeance was found to be consistent [START_REF] Falk | Subcontinuum mass transport of condensed hydrocarbons in nanoporous media[END_REF] with the so-called collective diffusivity 𝐷 0 = 𝐾 𝜌𝑘 B 𝑇 where 𝐷 0 is defined as Onsager's transport coefficient associated with the chemical potential gradient: 𝐽 = 𝜌𝑣 = -[𝜌𝐷 0 /(𝑘 B 𝑇)] ∇𝜇. This expression can be recovered by using Gibbs-Duhem equation at constant temperature, i.e. 𝜌d𝜇 = d𝑃.

The analogy between the permeance and collective diffusivity -albeit established for an in-16 O'Sullivan, June 1, 2022 compressible fluid -is an important cornerstone as it provides a means to estimate the permeance of a medium without performing sometimes complex permeability measurements. Of particular relevance, the collective diffusivity can be assessed using neutron scattering techniques with deuterated liquids (in contrast to hydrogen, deuterium is a strong neutron coherent scatterer which allows probing the collective part of the dynamical structure factor). A key advantage of such experiments over permeability measurements is that they are performed under no flow conditions -i.e. fluid at rest -as the fluctuation-dissipation theorem allows linking the transport coefficient (collective diffusivity here) to fluctuations in a given microscopic quantity (center of mass velocity here) through the Green-Kubo formalism. This central framework in statistical mechanics allows one to relate any transport quantity to fluctuations in microscopic variables [START_REF] Mc Quarrie | Statistical Mechanics[END_REF]. Of particular interest for fluid transport in confined media, it can be shown [START_REF] Barrat | Basic Concepts for Simple and Complex Liquids[END_REF] that the self-diffusivity is directly given by the time integral of the time autocorrelation function of the individual velocities (i.e. ⟨𝑣(𝑡)𝑣(0)⟩) while the collective diffusivity is given by the integral of the time autocorrelation of the center-of-mass velocities (i.e. ⟨𝑣 COM (𝑡)𝑣 COM (0)⟩).

UPSCALING PERMEABILITY

Despite decades of intense research by engineers and physicists, rigorous structure-property relations that link permeance (or hydraulic conductivity) and variables (fluid and material properties, thermodynamic conditions) for a given system are still lacking. In practice, the development of accurate and practical relationships has been hampered by the fact that the different communities interested in these issues often work independently from each other. On the one hand, researchers working on complex porous media such as in engineering fields mostly use Darcy's equation in an empirical fashion while trying to use reverse-engineering approaches to relate the microstructure of the porous medium to the flow properties observed in its pore space. On the other hand, researchers working in more basic science fields usually employ simpler porous media (ideally close to simple geometry channels) in combination with microscopic formalisms from statistical physics to unravel the complexity of fluid flow in confinement or near surfaces at the microscopic scale. 17 O'Sullivan, June 1, 2022

Engineering approaches: upscaling from pore scale

Full-field methods

Geomaterials contain solid grains and pore spaces of heterogeneous shapes, sizes and orientations. In full-field methods, features of the microstructure such as grains and pores are discretized in space, and field variables such as displacement, stress, pore pressure and porosity change are calculated at each point of the discretization grid. Composite materials were successfully simulated with the Finite Element Method (FEM), by discretizing each feature of the microstructure with several volume elements, so as to reproduce the spatial distribution of the composite hydro-mechanical properties in a statistical sense [START_REF] Chen | Numerical estimation of rev and permeability tensor for fractured rock masses by composite element method[END_REF]. The pore space can be modeled with cohesive zone elements (as a network of fractures) or with volume elements (as a porous solid with a shear modulus close to zero and with a permeability several orders of magnitude higher than that of the solid matrix around it). Alternatively, the FEM model can be designed to simulate fluid flow through a composite made of porous subdomains with different effective permeabilities as opposed to through the pore space itself. Displacements and pore pressures are solved at the nodes of the elements while stresses and porosity changes are calculated at integration points. If the nodes at the interface between subdomains are tied, fluid mass exchanges are calculated by assuming instantaneous equilibrium of pore pressures at the interface, i.e., infinite permeability orthogonal to the interface. Alternatively, interface elements can be inserted between the subdomains to reflect specific interface conductivity properties, including impermeable conditions. In all those FEM approaches, the field variables are not assumed to be uniform in each subdomain (or microstructure feature), which is why we call them "full-field methods". Their main limitation is that mesh generation remains a challenge for complex microstructure topologies. As noted above, all numerical methods such as FVM and LBM assume a given flow regime and boundary condition set whose valididty is clearly questionable in many instances. Fast Fourier Transform (FFT) -based methods are attractive for highly heterogeneous materials, because calculations are carried out on a regular grid and periodic boundary conditions minimize size effects. The numerical method established by [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] is based on the FFT and on the introduction of a reference mate-18 O'Sullivan, June 1, 2022 rial leading to the Lippmann-Schwinger (LS) equation [START_REF] Suquet | Une méthode simplifiée pour le calcul des propriétés élastiques de matériaux hétérogènes à structure périodique[END_REF], which is then discretized.

Variational estimates of poromechanical and permeability properties were derived from the LS equation [START_REF] Bignonnet | Fft-based homogenization of permeability using a hashin-shtrikman type variational framework[END_REF][START_REF] Brisard | Variational estimates of the poroelastic coefficients[END_REF]. Full-field approaches based on the LS equation typically perform better than a FEM [START_REF] Lebensohn | Orientation imagebased micromechanical modelling of subgrain texture evolution in polycrystalline copper[END_REF], but by construction, FFT-based methods are better suited when the microstructure can be represented by periodic cells. For example, the permeability of a periodic micro-porous medium saturated by a viscous fluid was successfully calculated by using a FFT based computing method in which slip boundary conditions at the pore wall are accounted for by introducing an interphase of inifinitely small thickness between the fluid and the solid [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF]. Later, a variational FFTbased scheme based on the minimization of an energy Hashin-Shtrikman functional was proposed, which allowed calculation of rigorous bounds of homogenized elastic properties in an energetically consistent way [START_REF] Brisard | Fft-based methods for the mechanics of composites: A general variational framework[END_REF]. This variational FFT approach makes it possible to perform computations on large discretization grids [START_REF] Brisard | Combining galerkin approximation techniques with the principle of hashin and shtrikman to derive a new fft-based numerical method for the homogenization of composites[END_REF]. The Hashin and Shtrikman -based FFT scheme then inspired a numerical scheme to solve the Darcy flow problem with periodic boundary conditions and calculate the rigorous bounds for the permeability tensor [START_REF] Bignonnet | Fft-based bounds on the permeability of complex microstructures[END_REF]. One of the key issues for random composites is to find the size of a statistically representative volume over which to calculate the effective properties. An automated procedure to define the size of the discretized domain was recently proposed by [START_REF] Bignonnet | Efficient fft-based upscaling of the permeability of porous media discretized on uniform grids with estimation of rve size[END_REF], who found that the discretization of trial force fields on a regular grid permits the homogenization of permeability from Stokes equations in voxel representations of the pore space and that energy consistent discretizations deliver the most accurate results. Of note, the FFT based computing method proposed in [START_REF] Monchiet | A fft-based method to compute the permeability induced by a stokes slip flow through a porous medium[END_REF] to solve the Stokes equations was extended to doubly porous solids constituted of two populations of pores which are different in size [START_REF] Mezhoud | Computation of macroscopic permeability of doubly porous media with fft based numerical homogenization method[END_REF].

Pore network models

Pore network models (PNMs) (Figure 1 (d)) can be used to relate the permeability tensor to microstructure geometric properties. The original Kozeny-Carman model [START_REF] Kozeny | Uber kapillare leitung der wasser in boden[END_REF][START_REF] Carman | Fluid flow through granular beds[END_REF] considers that the pore fluid flows in a bundle of parallel pipes. Modified Kozeny-Carman permeability formulas [START_REF] Brace | Permeability of granite under high pressure[END_REF][START_REF] Berryman | Use of digital image analysis to estimate fluid permeability of porous materials: Application of two-point correlation functions[END_REF][START_REF] Petersen | Relations between specific surface area and soil physical and chemical properties[END_REF][START_REF] Mavko | The effect of a percolation threshold in the kozeny-carman relation[END_REF][START_REF] Saar | Permeability-porosity relationship in vesicular basalts[END_REF] indirectly relate average flow properties to tortuosity (e.g., through electrical conductivity). Statistical flow networks models [START_REF] Dienes | Permeability, percolation and statistical crack mechanics[END_REF][START_REF] Gueguen | Transport properties of rocks from statistics and percolation[END_REF][START_REF] Torquato | Relationship between permeability and diffusion-controlled trapping constant of porous media[END_REF][START_REF] Schubnel | Quantifying damage, saturation and anisotropy in cracked rocks by inverting elastic wave velocities[END_REF][START_REF] Fortin | Influence of thermal and mechanical cracks on permeability and elastic wave velocities in a basalt from mt. etna volcano subjected to elevated pressure[END_REF][START_REF] Arson | Influence of damage on pore size distribution and permeability of rocks[END_REF] and fractal network models [START_REF] Tyler | Fractal processes in soil water retention[END_REF][START_REF] Neuman | On advective transport in fractal permeability and velocity fields[END_REF][START_REF] Pape | Permeability prediction based on fractal pore-space geometry[END_REF] depend on statistical microstructure descriptors. In the percolation theory [START_REF] Stauffer | Introduction to percolation theory[END_REF], the characteristic flow path length is postulated, and transitions in flow regimes can only be predicted if the probability of fluid saturation is known for each pore or crack. Although a non-uniform pore arrangement produces hysteresis in retention curves [START_REF] Romero | Water permeability, water retention and microstructure of unsaturated compacted boom clay[END_REF][START_REF] Wheeler | Coupling of hydraulic hysteresis and stressstrain behaviour in unsaturated soils[END_REF][START_REF] Nuth | Advances in modelling hysteretic water retention curve in deformable soils[END_REF], most models that relate capillary pressure to pore size [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF][START_REF] Pereira | Retention and permeability properties of damaged porous rocks[END_REF] assume that the pore network is a bundle of pipes of constant cross section, which are entirely filled with either wetting or non-wetting fluid.

More recent PNMs comprise pores (nodes) connected by constrictions or throats (edges).

A partitioning algorithm is required to identify individual pores if the pore space topology is determined from micro computed tomography [START_REF] Dong | Pore-network extraction from micro-computerized-tomography images[END_REF][START_REF] Taylor | A new method to identify void constrictions in micro-ct images of sand[END_REF] or from a DEM model [START_REF] Reboul | A computational procedure to assess the distribution of constriction sizes for an assembly of spheres[END_REF][START_REF] Van Der Linden | A computational geometry approach to pore network construction for granular packings[END_REF]. Partitioning is non-trivial as the void space is continuous and no partitioning algorithm can be wholly objective. For example, [START_REF] Sufian | Ability of a pore network model to predict fluid flow and drag in saturated granular materials[END_REF] compare PNMs obtained using slightly different approaches to partitioning.

The model then considers mass balance for each pore in the system:

Σ 𝑏 𝑁𝑝 𝑏=𝑏 1 𝑄 𝑎𝑏 = 0 (12) 20 O'Sullivan, June 1, 2022
where 𝑄 𝑎𝑏 is the fluid flow rate from pore 𝑎 to pore 𝑏 and 𝑏 𝑗 , 𝑗 = 1 -𝑁 𝑝 are the ids of the 𝑁 𝑝 pores adjacent to pore 𝑎. Darcy's law effectively governs the flow from one pore to the next as

𝑄 𝑎𝑏 = 𝐺 𝑎𝑏 𝑝 𝑎 -𝑝 𝑏 𝑙 𝑎𝑏 (13)
where 𝐺 𝑎𝑏 is the hydraulic conductance of the throat connecting pores 𝑎 and 𝑏; the intrinsic conductance is given by 𝐾 𝑎𝑏 = 𝜇𝐺 𝑎𝑏 . The proposed models to determine 𝐾 𝑎𝑏 use a measure of the pore throat geometry and consider the Hagen-Poiseille equation; expressions have been proposed

by [START_REF] Bryant | Prediction of relative permeability in simple porous media[END_REF], [START_REF] Chareyre | Pore-scale modeling of viscous flow and induced forces in dense sphere packings[END_REF] and [START_REF] Patzek | Shape factor and hydraulic conductance in noncircular capillaries: I. one-phase creeping flow[END_REF].

Mean-field homogenization methods

In the homogenization theory, the Representative Elementary Volume (REV) is defined as the volume that has a characteristic size equal to the length of the shortest path that needs to be traveled inside the continuum to observe uniform variability of permeability properties. Typically, the REV is two orders of magnitude larger than the features of the microstructure that are expected to influence the permeability of the medium (e.g., pores, cracks). In mean-field homogenization approaches, the REV is divided into subdomains that each represent a microstructure feature (inclusion or matrix in which the inclusions are embedded in), and the field variables such as stress and pore pressure are uniform in each subdomain. This is in contrast to full-field methods, in which each feature or subdomain is discretized in space, and to PNMs, in which the concept of a matrix does not exist and the local gradient of pore pressure needs to be calculated at every constriction. Most permeability models based on mean-field homogenization approaches consider that Hagen-Poiseuille flow takes place in inclusions that represent penny-shaped cracks that are connected to a porous matrix and can intersect with one another. The matrix/inclusion interaction equation for mass transport is [START_REF] Dormieux | Approche micromécanique du couplage perméabilitéendommagement[END_REF]:

∇𝑝(x) = A(x) . ⟨∇𝑝(x)⟩ 𝐷 = A(x) . ∇ P (14) 21 O'Sullivan, June 1, 2022
in which ⟨...⟩ 𝐷 is the space average over the domain 𝐷 of the REV, 𝑝(x) and P are the local and REV scale pore pressure fields, and A(x) is the localization tensor (see Figure 3). Since the average velocity in the REV is equal to the average of the local velocity fields in the matrix and inclusions:

k . ∇ P = ⟨k(x).∇𝑝(x)⟩ 𝐷 (15)
where k and k(x) are the macroscopic (REV scale) and microscopic permeability tensors, respectively. The knowledge of A(x) yields the expression of the macroscopic permeability tensor:

k = ⟨k(x).A(x)⟩ 𝐷 (16)
The expression of A(x) was found for Eshelby's self-consistent scheme [START_REF] Dormieux | Approche micromécanique du couplage perméabilitéendommagement[END_REF][START_REF] Barthélémy | Effective permeability of media with a dense network of long and micro fractures[END_REF][START_REF] Chen | Micromechanical modeling of anisotropic damageinduced permeability variation in crystalline rocks[END_REF] and for Ponte-Castenada-Willis variational approach [START_REF] Levasseur | A micro-macro approach of permeability evolution in rocks excavation damaged zones[END_REF], which was used to predict the evolution of damage and permeability in the excavation damaged zone. More recently, [START_REF] Abdalrahman | Is trabecular bone permeability governed by molecular ordering-induced fluid viscosity gain? arguments from re-evaluation of experimental data in the framework of homogenization theory[END_REF] used a self-consistent approach to express the localization tensor of interpenetrating cylindrical pore channels oriented randomly and interacting with impermeable spherical inclusions. The model was applied to the prediction of bone permeability, which is critical to fracture healing. The self-consistent model was then used to understand fluid flow -driven bone stimulation at three scales [START_REF] Estermann | Multiscale modeling provides differentiated insights to fluid flow-driven stimulation of bone cellular activities[END_REF].

Permeability upscaling assisted by Artificial Intelligence (AI)

In geomechanics, feature selection algorithms have been employed to identify sets of pore network and microstructural properties that optimally characterize permeability [START_REF] Van Der Linden | Machine learning framework for analysis of transport through complex networks in porous, granular media: A focus on permeability[END_REF]. Long-Short-Term-Memory (LSTM) neural networks have been used to optimize the expressions of the meso-and macro-effective stress and permeability tensors from offline simulations at the micro-and meso-scales, respectively [START_REF] Wang | A multiscale multi-permeability poroplasticity model linked by recursive homogenizations and deep learning[END_REF]. Hydraulic soil responses can also be generated from a neural network trained "off-line" with LBM results: mechanical constitutive updates are obtained from DEM simulations linked to a macroscopic dual-permeability model 22 O'Sullivan, June 1, 2022 [START_REF] Wang | An updated lagrangian lbm-dem-fem coupling model for dual-permeability fissured porous media with embedded discontinuities[END_REF]. Hybrid algorithms can predict either or both the mechanical and hydraulic behavior of poroelastic media via constitutive models or via Machine Learning (ML) [START_REF] Bahmani | A kd-tree-accelerated hybrid data-driven/model-based approach for poroelasticity problems with multi-fidelity multi-physics data[END_REF]. The ML algorithm solves an optimization problem (finding state variables as close as possible to the data set) under constraints (momentum/mass balance equations, boundary conditions). Reinforced learning allowed optimizing relationships between microstructure descriptors and state variables (Wang and Sun, 2019a;Wang et al., 2019). The problem effectively is set as a board game, in which a finite number of microstructure descriptors and state variables form the nodes of a directed graph. The edges represent yet-unknown definitions or constitutive laws. The probability of choosing an edge versus another is governed by a reward system. The reward is calculated based on the accuracy of the constitutive model under constraints applied on the form of the graph. Self-improving algorithms learn not only the constitutive laws, but also the optimal setup, architecture, and hyper-parameters of the deep Convolutional Neural Networks (CNNs). The deep reinforcement learning strategy is based on a Markov-Decision process and was applied to predict path-dependent retention curves and anisotropic flow in micropores [START_REF] Heider | An offline multi-scale unsaturated poromechanics model enabled by self-designed/self-improved neural networks[END_REF].

Statistical physics approaches: from molecular to pore scale

In order to bridge the gap between molecular dynamics and transport in porous networks, it is necessary to: (i) characterize the structure of the porous medium and measure transport properties at multiple scales, and (ii) couple experimental and numerical techniques to determine the expression of the pore network permeability (or permeance) as a function of geometric and physical properties established at multiple scales. In this context, statistical physics approaches -which were briefly mentioned and introduced in this paper -are necessary to go beyond the classical Stokes equation and include additional effects that induce deviations from Darcy's law. Such a powerful framework is also needed to provide a more general Darcy's scheme, where the (intrinsic) permeability is replaced with complex expressions for the general permeance.

Geometric and transport properties at different scales

To provide a robust scheme for transport upscaling, experimental structural data must be gathered at each relevant scale. Such characterization methods, which are listed here in a non-exhaustive 23 O'Sullivan, June 1, 2022 fashion, allow probing additional effects that may be taken into account by operating at different length scales. This includes techniques such as:

• Neutron or X-ray scattering at the ∼ nm scale and beyond to probe adsorption and, in general, all physical and physicochemical effects;

• Transmission Electron Microscopy (TEM), electron energy loss spectroscopy (EELS), Nanotomography (NanoT-SCAN) at the ∼ 10 nm scale to probe pore shape effects;

• Focused ion beam and regular scanning electron microscopy (FIB-SEM and SEM) to probe pore networks effects at the ∼ 100 nm -1 𝜇m scales;

• Computed tomography (CT-SCAN) and conventional microscopy at the macroscopic scales.

Coupling experiments and simulations to find the expression of permeability

Experimentally, many techniques allow probing the dynamics of fluids at different scales, e.g., relaxometry nuclear magnetic resonance (NMR), pulsed field gradient nuclear magnetic resonance (PFGNMR), quasielastic neutron scattering (QENS), permeability measurements on core samples.

Similarly, several numerical methods can be used, such as molecular dynamics simulation (MD), brownian dynamics (BD), dissipative particle dynamics (DPD), smoothed particle hydrodynamics (SPH), lattice Boltzmann method (LBM), finite difference (FD) or finite volume (FV) reservoir simulators. Connections between these different techniques could pave the way for robust upscaling, whereby details of the structure-flow relationship are upscaled at the larger scale while keeping as much information as possible. A key step along this roadmap is to collect transport data using other techniques than permeability experiments. Indeed, this is needed as permeability (or permeance)

is the key quantity that should be modeled/predicted in the end so that it cannot be used directly as a variable in the upscaling approach. For example, the collective diffusivity 𝐷 0 -which can be assessed using techniques such as neutron scattering -is directly linked to the permeance 𝐾

PERSPECTIVES

In engineering practice Darcy's law is most often applied empirically to characterize transport in complex porous media, where permeability is more a descriptor (a tensor when non-isotropic materials are considered such as in sedimentary rocks and natural soils) rather than an intrinsic physical constant of the medium. Even when a purely macro-scale, empirical approach is adopted, it is challenging to accurately measure permeability in situ or in the laboratory [START_REF] Cashman | Groundwater Lowering in Construction A Practical Guide to Dewatering[END_REF]. Simple expressions to determine the permeability k based on flow and material parameters are currently unavailable for real, complex media (despite attempts to use physically robust approaches such as mean-field methods).

This paper has highlighted that recent developments in numerical methods, imaging and advances in computational power enable us to take a new perspective on flow in porous media and advance understanding of the fundamentals that underlie Darcy's law. However, the paper has also shown the benefit of using an analytical approach to explore the limitations of the applicability of Darcy's law in a very direct way. Starting from a consideration of flow in a single, cylindrical pore, it has explored the applicability of Darcy's law in different scenarios, including non-laminar flow and non-advective flow. Working from this idealized model has enabled a reflection on the possible definitions of permeability to extend Darcy's law to poromechanical problems that involve strong couplings with chemical and physical processes that occur at the pore scale and below. Using this model, it can be shown that Darcy's law is valid at all scales where the fluid can be modelled as a continuum so long as the permeability is expressed as the correct function of geometric, chemical, physical and mechanical parameters that fluid flow depends on. The exact dimensions at which the continuum assumption applies depends on the fluid in question, the temperature and the pressure; in general, for fluids, the continuum-discrete limit is approximately between 1nm and a few nm 25 O'Sullivan, June 1, 2022 [START_REF] Bocquet | Nanofluidics, from bulk to interfaces[END_REF]. At the continuum scale, the value of considering a model of single cylindrical pore, albeit a highly simplified abstraction of the situation in a real material, is that it provides us with a clear basis to understand some of the limitations of the applicability of Darcy's law.

In some situations though, expressing permeability as an appropriate function remains a challenge. To date, there are still many unexpected instances in which fluid velocity is not a linear function of the gradient of pore pressure. From a geotechnical engineering perspective, for example, there are some questions on the validity of Darcy's law in clays at low pressure gradients because the structure of water in very small pores may be an issue [START_REF] Hansbo | Consolidation equation valid for both darcian and non-darcian flow[END_REF]. This is an issue that must be addressed using simulations at the molecular scale (i.e. Molecular Dynamics) and it is timely to harness this type of numerical model to look at this issue. We also need to better understand how Darcy's law can be adapted to predict the infiltration of suspensions of bentonite/polymers used to support boreholes and excavations: in these cases flow depends on the rheology of the fluid, which is non-Newtonian [START_REF] Ejezie | Permeation behaviour of phpa polymer fluids in sand[END_REF]; Darson-Balleur ( 2019)).

One of the key challenges is to characterize microstructure and measure transport properties at multiple scales. To measure transport coefficients for different gases at various temperatures and densities and establish accurate structure-property relationships, the impact of thermodynamic and transport conditions onto fluid flow must be better understood. In particular, despite their broadly acknowledged contribution to fluid flow, adsorption phenomena but also other thermodynamical mechanisms (e.g. condensation, cavitation, pore blocking, etc.) are often overlooked as they are treated in an effective, coarse-grained fashion. In this context, while these effects are often assumed to occur only in very confining materials (e.g. nanoporous solids), they can have a strong impact even when large pore materials are considered. For instance, even a vanishing amount of such nanoporosity in the host medium can induce drastic macroscopic effects by changing locally the thermodynamic and hydrodynamic boundary conditions imposed on the large porosity. Establishing a robust connection between permeance (or hydraulic conductivity) and the variables of the problem for simple pores is a prerequisite to any upscaling attempt to describe transport in complex, multi- (1) Averaging. In the case of a domain 𝐷 that is subjected to a linear pore pressure in the far field (𝑃 ∞ (x) = x • ∇𝑃 ∞ ) and that contains a heterogeneity Ω, the gradient of pressure in the matrix 𝐷\Ω, noted ∇𝑝 0 (x), is equal to the uniform gradient of pressure applied at the far-field boundary (∇𝑃 ∞ ) plus the gradient of a disturbance pressure field (∇𝑝 𝑑 (x)), due to the presence of the heterogeneity. Similarly, the gradient of pressure in the heterogeneity (∇𝑝 𝑖 (x)) is equal to ∇𝑃 ∞ + ∇𝑝 𝑑 (x). Additionally, the flow rates in the matrix (v 0 (x)) and in the heterogeneity (v i (x)) are related to 𝑝 0 (x) and 𝑝 𝑖 (x) respectively, via Darcy's law. In Eshelby's theory, the gradients of the local fields of pressure are related to the average gradient of pressure in the domain 𝐷 by localization tensors (A 0 (x) for the matrix, A i (x) for the heterogeneity), which depend on the shape of the heterogeneities.

(2) Solving for local field variables. The pressure field in the matrix, 𝑝 0 (x), is calculated by solving a cavity expansion problem, in which 𝑃 ∞ (x) is applied in the far field and 𝑝 𝑖 (x) is applied at the matrix/heterogeneity interface. The field 𝑝 𝑖 (x) is calculated by the equivalent inclusion method, in which the heterogeneity Ω is assigned the same permeability as the matrix 𝐷\Ω. To compensate the error made in the estimation of the flow rate in the heterogeneity, a gradient of eigenpressure (∇𝑝 * 𝑖 (x)) is subtracted from the local gradient of pressure (∇𝑝 𝑖 (x)) in Darcy's law. It can be shown that the gradient of eigenpressure field is uniform in Ω and that there is a linear relationship between the gradient of the disturbance pressure field (∇𝑝 𝑑 (x)) and the gradient of the eigenpressure field (∇𝑝 * 𝑖 ). The relationship between the two holds in the so-called Eshelby tensor S i , which itself, depends on the Green's tensor and on the localization tensors (A 0 (x), A i (x)). (3) Generalization. For a set of heterogeneities, assumptions on the interactions between heterogeneities are introduced to close the formulation.

•

  A(x) localization tensor • 𝑏 slip length • 𝑏 0 system dependant constant • 𝐶 constant of integration • 𝑐 sound velocity (scalar) in fluid • 𝐶𝑁 𝑁s Convolutional Neural Networks • 𝐷 diameter • 𝐷 0 collective diffusivity • DPD dissipative particle dynamics • FVM Finite Volume Method • FD finite difference • g acceleration due to gravity • 𝐺 𝑎𝑏 hydraulic conductance of throat connecting pores 𝑎 and 𝑏 • J flux • 𝑗 𝑥 (k, 𝑡) Fourier component of the flow velocity • k permeability tensor • k(x) local permeability tensor • K permeance tensor • K ′ hydraulic conductivity • 𝐾 𝑎𝑏 intrinsic conductance • 𝑘 𝐵 Boltzmann's constant • Kn Knudsen number • 𝐿 measure of length • LBM Lattice Boltzman Method • LSTM Long-Short-Term-Memory • Ma Mach number 2 O'Sullivan, June 1, 2022

•

  𝑡 time • v velocity field • v average flow rate • v s seepage velocity • 𝜂 fluid viscosity • 𝜆 mean free path length • 𝜉 friction between the flowing fluid and solid • 𝜌 fluid density • 𝜎 kinetic diameter • 𝜎 𝛼𝛽 stress tensor component • 𝜏 𝑅 stress relaxation time in the fluid particle • 𝜏 fluid momentum transfer time INTRODUCTION In his pioneering work on city fountains in Dijon (France), Henry Darcy established an experimental relationship in which the average flow rate 𝑣 is proportional to the pressure drop, Δ𝑃,

  our understanding of flow at different scales within the framework of Darcy's law. Three questions underpin this discussion: (1) How can we establish the validity of Darcy's law at multiple scales? (2) What mechanisms underlie the regimes where Darcy's law fails to capture the relationship between the flow rate and the pressure gradient? (3) How can we upscale permeability from fundamental measures of porous media? This contribution adopts a fundamental perspective. It does not address key issues in engineering practice such as challenges associated with measuring permeability in situ or in the laboratory. However, the insight emerging from the detailed consideration of the nature of fluid flow in porous media can underpin practical analyses and interpretation of in situ observations.

  , Zhao and O'Sullivan (2022)), laser scanning of transparent soil (e.g. Sanvitale et al. (2022)), or micro computed tomography scanning (e.g. Talon et al. (2012); Blunt et al. (2013); Taylor et al. (2015); Sun and Wong (2018)). The particle surfaces are taken as boundaries to the fluid flow. In this type of model, rather than considering Darcy's law, the Navier-Stokes equations are solved numerically using, for example, the finite volume method (FVM) (Figure 1 (b)). The Lattice-Boltzmann Method (LBM) which solves Boltzmann's transport equation by introducing a lattice discretization can also be used (Figure 1 (c))

  Re = 𝜌𝑣𝐷/𝜂 with 𝜌 the fluid density, 𝑣 the flow rate, 𝐷 the characteristic pore size and 𝜂 the fluid dynamic viscosity). In the case of larger Reynolds numbers, this approximation no longer holds and 11 O'Sullivan, June 1, 2022 quadratic terms in the flow rate must be added to Darcy's equation. In what follows, we consider

  , June 1, 2022 where 𝜉 is the friction between the flowing fluid and the solid. Writing this conservation relation in Fourier space in a direction 𝑥 normal to the flow direction gives the following differential equation 𝜕 𝑗 𝑥 (k, 𝑡)/𝜕𝑡 = (-𝜂𝑘 2 /𝜌 -𝜉) 𝑗 𝑥 (k, 𝑡) where 𝑗 𝑥 (k, 𝑡) is the Fourier component of the flow velocity.

  [START_REF] Bocquet | Nanofluidics, from bulk to interfaces[END_REF] provide an explanation for the viscous flow approximation breakdown, which is summarized here. Continuum conservation laws such as Stokes' equation assume a time scale separation between the stress relaxation time 𝜏 𝑅 inside the fluid particle and the fluid momentum transfer time 𝜏. On the one hand, stress relaxation within the fluid occurs over a typical time given by the time correlations in the stress tensor since 𝜂 = 1/𝑉 𝑘 B 𝑇 ∫ 𝜎 𝛼𝛽 (𝑡)𝜎 𝛼𝛽 (0)d𝑡; in this equation, 𝜎 𝛼𝛽 is the stress tensor component with 𝛼, 𝛽 = 𝑥, 𝑦, or 𝑧. This expression, known

  The Darcy/Stokes coupled equations are solved by introducing the Brinkman equation to represent the fluid exchange at the interface between the different subdomains. The fields of velocity and traction are continuous across the interface, which is compatible with the use of the Fourier series. Despite these advances, the reconstruction of a kinematically admissible velocity 19 O'Sullivan, June 1, 2022 field from the force field and the coupling between fluid flow and pore deformation remain open issues, even in the most recent FFT-based methods.

[

  𝐾 ∼ 𝐷 0 /𝜌𝑘 B 𝑇] as noted earlier. Finally, we note that most upscaling techniques (e.g. pore network models) and simulation techniques (e.g. LBM) -despite being established on robust grounds -only include in an effective fashion the complex thermodynamics of the confined fluids 24 O'Sullivan, June 1, 2022 (adsorption, condensation, cavitation, etc). In contrast, other lattice models (e.g. Boţan et al. (2015)) which account for the change in the chemical potential of the fluid upon stationary and/or transient transport regimes allow one (in principle) to describe the coupling between fluid transport and thermodynamics in complex porous media.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Sub-pore scale modeling of fluid flow. (a) Poiseuille flow for a fluid made of molecules (red spheres) in a slit pore of width 𝐷 (stone like aspect). The fluid flow obeys a parabolic velocity profile upon applying a pressure gradient ∇𝑃 parallel to the pore surfaces. 𝑥 and 𝑟 are the direction parallel and perpendicular to the pore surface. (b) Depending on the so-called surface boundary condition, 𝑣 𝑥 (𝑟) ∼ 0 at the pore surface (no-slip boundary condition, left panel) or 𝑣 𝑥 (𝑟) ≠ 0 (slip boundary condition, right panel). In the latter case, slippage occurs with a slip length 𝑏 which corresponds to the distance at which the extrapolated velocity profile reaches zero: -𝑏𝜕𝑣 𝑥 /𝜕𝑟 = 𝑣 𝑥 .

  ∼ ∇𝑃 and that the permeability 𝑘 depends on the fluid molecule as well as on the fluid density. The concept of permeability can encompass both geometrical effects (tortuosity, morphology), physical/chemical effects (e.g. adsorption) and mechanical effects (non viscous effects, necking, etc.)

[START_REF] Falk | Subcontinuum mass transport of condensed hydrocarbons in nanoporous media[END_REF]

who considered hydrocarbon transport (e.g. methane, propane, hexane, nonane, and dodecane) in a disordered nanoporous solid. Falk et al. found that the flow varies linearly 15 O'Sullivan, June 1, 2022 with the prescribed pressure gradient i.e. J

TABLE 1 .

 1 Typical values for hydraulic conductivity (also called coefficient of permeability in the geotechnical literature) for geomaterials developed with reference to[START_REF] Cashman | Groundwater Lowering in Construction A Practical Guide to Dewatering[END_REF]

	Soil Type	Hydraulic	conductivity	Qualitative / engineering de-
		K ′ = k𝜌𝑔/𝜂 (𝐾 ′ 𝑖 𝑗 in 𝑚𝑠 -1 )	scription
	Intact clays	< 10 -9		Practically impermeable
	Fissured or laminated clays	10 -9 to 10 -7		Very low
	Sandy silts, very silty fine	10 -8 to 10 -5		Low to very low
	sands, and laminated or mixed
	strata of silt/sand/clay
	Silty sands	10 -6 to 10 -4		Low
	Fine and medium sands	10 -4 to 5 × 10 -4		Moderate to low
	Clean sand and sand/gravel	5 × 10 -4 to 10 -3		High to moderate
	mixtures
	Clean gravels	> 10 -3		High

O'Sullivan, June 1, 2022

ACKNOWLEDGEMENTS

Some of the arguments made in this paper are inspired from a formal debate held by COS and BC and moderated by CA during the Biot-Bažant Conference on Engineering Mechanics and Physics of Porous Materials, held June 1-3, 2021. A recording of the debate is available in (O' Sullivan and Coasne, 2021). We thank R. Hartkamp for providing the data for the molecular simulation snapshot in Figure 1 (a). We also thank Mr. Tokio Morimoto who provided the images in Figures 1 (b) and (d). Some of the arguments presented in this paper originated from Chloé Arson's project "BRITE Pivot Track: Micro-macro modeling of reactive flow and rock weathering enhanced by Artificial Intelligence" which is funded by the U.S. National Science Foundation (NSF) under grant CMMI # 2135584. The financial support from the NSF is gratefully acknowledged.

O'Sullivan, June 1, 2022 scale systems. Software development is still needed to decipher the scales of interest in fluid flow in porous media. Key challenges remain, however there is now a real possibility of addressing some complex fundamental questions using recent developments in numerical algorithms and advances in computational power. This is an area which can be explored by appropriate application of direct numerical simulation using either FVM or LBM and potentially including pore network modelling.

DATA AVAILABILITY STATEMENT

No data, models, or code were generated or used during the study. (c) Flow velocity field in this pore space, e.g. calculated by LBM (Lattice Boltzman Method) after [START_REF] Su | Exploring semi-solid alloy deformation with discrete element method simulations and synchrotron radiography[END_REF]. In the case of a domain 𝐷 that is subjected to a linear pore pressure in the far field (𝑃 ∞ (x) = x • ∇𝑃 ∞ ) and that contains a heterogeneity Ω, the gradient of pressure in the matrix 𝐷\Ω, noted ∇𝑝 0 (x), is equal to the uniform gradient of pressure applied at the far-field boundary (∇𝑃 ∞ ) plus the gradient of a disturbance pressure field (∇𝑝 𝑑 (x)), due to the presence of the heterogeneity. Similarly, the gradient of pressure in the heterogeneity (∇𝑝 𝑖 (x)) is equal to ∇𝑃 ∞ + ∇𝑝 𝑑 (x).

Additionally, the flow rates in the matrix (v 0 (x)) and in the heterogeneity (v i (x))

are related to 𝑝 0 (x) and 𝑝 𝑖 (x) respectively, via Darcy's law. In Eshelby's theory, the gradients of the local fields of pressure are related to the average gradient of pressure in the domain 𝐷 by localization tensors (A 0 (x) for the matrix, A i (x) for the heterogeneity), which depend on the shape of the heterogeneities.

(2) Solving for local field variables. The pressure field in the matrix, 𝑝 0 (x), is calculated by solving a cavity expansion problem, in which 𝑃 ∞ (x) is applied in the far field and 𝑝 𝑖 (x) is applied at the matrix/heterogeneity interface. The field 𝑝 𝑖 (x) is calculated by the equivalent inclusion method, in which the heterogeneity Ω is assigned the same permeability as the matrix 𝐷\Ω. To compensate the error made in the estimation of the flow rate in the heterogeneity, a gradient of eigenpressure (∇𝑝 * 𝑖 (x))

is subtracted from the local gradient of pressure (∇𝑝 𝑖 (x)) in Darcy's law. It can be shown that the gradient of eigenpressure field is uniform in Ω and that there is a linear relationship between the gradient of the disturbance pressure field (∇𝑝 𝑑 (x))

and the gradient of the eigenpressure field (∇𝑝 * 𝑖 ). The relationship between the two holds in the so-called Eshelby tensor S i , which itself, depends on the Green's tensor and on the localization tensors (A 0 (x), A i (x)).