Local Thermodynamic Description of Isothermal Single-Phase Flow in Simple Porous Media

Darcy’s law for porous media transport is given a new local thermodynamic basis in terms of the grand potential of confined fluids. The local effective pressure gradient is determined using non-equilibrium molecular dynamics, and the hydraulic conductivity and permeability are investigated. The transport coefficients are determined for single-phase flow in face-centered cubic lattices of solid spheres. The porosity changed from that in the closest packing of spheres to near unity in a pure fluid, while the fluid mass density varied from that of a dilute gas to a dense liquid. The permeability varied between 5.7×10-20m2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5.7 \times {10^{-20}} \hbox {m}^2$$\end{document} and 5.5×10-17m2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5.5 \times {10^{-17}} \hbox {m}^2$$\end{document}, showing a porosity-dependent Klinkenberg effect. Both transport coefficients depended on the average fluid mass density and porosity but in different ways. These results set the stage for a non-equilibrium thermodynamic investigation of coupled transport of multi-phase fluids in complex media.


Introduction
Porous media are everywhere in nature and technology, and transport through them is important. To take some widely different examples; we need to describe the transport of nanoparticles across cell layers with medicine to cancerous tissue (Stylianopoulos et al. 2018). We also need to describe the selective transport across the porous separators in batteries and fuel cells (Zlotorowicz et al. 2015). We are aiming for a description that reflects the underlying properties of the single pores.
Recent developments in imaging techniques and computer capabilities have improved our understanding of the physics at pore-scale considerably (Blunt 2017). But there is no consensus on how to upscale from the pore-scale to the Darcy scale. On the Darcy scale, transport is described as taking place between representative elementary volumes (REVs).
The REV in this work is defined to be large enough to be statistically representative of the system. From a statistical mechanics point of view, the REV includes all available microstates of the system. In this work, we are investigating a compressible single-phase fluid in a porous structure made up of solid particles in a face-centered cubic (fcc) lattice. The system is illustrated in Fig. 1. The blue large spheres represent solid, while the red particles represent fluid. We shall use additive variables to define the REV, similar to Whitaker (1986). This definition of the REV differs from the one suggested by Nordahl and Ringrose (2008) who used a constant permeability as a criterion. Because of the fcc symmetry, a unit cell is here a proper choice of REV. Around any point in the porous medium, we can choose a unit cell with this point at its center. In this way, we can obtain a continuous path between unit cells on the REV scale. We first determine the pressures of the REV at equilibrium, as a function of the temperature, fluid mass density, and porosity, p(T, f , ) . This can be regarded as finding the equation of state of the REV. This is analogous to equations of state for the capillary pressure in two-phase flow, where the capillary pressure is, for example, assumed to be a function of the saturation of one component (Hassanizadeh and Gray 1993). Up to this point, we have used REV densities much like Whitaker (1986).

Fig. 1
The representative elementary volume (REV) for a single-phase fluid in a facecentered cubic lattice of solid spheres is the same size as the unit cell of the lattice. The radius of the solid particles is R and the lattice constant is a. The squares with dashed lines mark the magnitude of the REV. All regions have the same thermodynamic properties Thermodynamic theories for transport in porous media are little developed. In confinement, fluid thermodynamic properties will deviate from their corresponding bulk phase values. In an attempt to find a continuous description on the Darcy scale, we need to reflect on the properties of the pore scale or sub-pore scale, including the nanometer scale. Hence, we need descriptions of fluxes and forces for the REV of the porous medium. Is it possible to find a description that does not explode in complexity, but still brings forward the characteristic properties of the smaller scale, such as wetting and adsorption? Pore-scale descriptions of REVs have contained up to 27 variables (Gray and Hassanizadeh 1998). We have claimed that we can reduce the number of variables to a more practical number, and describe here the first steps in the direction to apply this relatively new thermodynamic procedure for coarse-graining Rauter et al. 2020;Bedeaux and Kjelstrup 2021).
In our search for a thermodynamic coarse-graining procedure to be applied to transport in porous media, we have chosen to define effective thermodynamic variables by sets that combine in an additive manner (Kjelstrup et al. 2018Bedeaux and Kjelstrup 2021). We have derived the Gibbs equation for the REV, and in the analysis of transport, we have assumed the validity of the Gibbs equation. This gives the entropy production for the REV and a thermodynamic basis for Darcy's law in an isothermal system. Our proposal to find the equation of state and the Gibbs equation is new. Here, we apply the method to isothermal, compressible, single-phase flow. This entails obtaining the REV grand potential as the sum of contributions from all phases and surfaces. A theory of transport on the Darcy scale follows in a way that is standard in non-equilibrium thermodynamics.
An excellent tool to analyze transport in porous media is non-equilibrium molecular dynamics (NEMD) (Todd and Daivis 2017). This tool allows us to simulate molecular properties (like velocities and forces), yet upscale to fluid properties (like pressure) on the Darcy scale. In NEMD, we solve Newton equations for particles, so the outcome can also be used to assess assumptions made in the thermodynamic theory. In this sense, the tool supplements lattice Boltzmann simulations and numerical solutions to the Navier-Stokes equation. A downside is that the length and time scale become limited as NEMD is a computationally expensive technique. Here, we will use NEMD to simulate the flow of a model system, with methane-like molecules in a face-centered cubic (fcc) lattice made up of spherical solid particles of varying porosity, see Fig. 2. We will investigate a vast range of fluid densities, which varies from highly compressible to nearly incompressible. NEMD has been used to simulate many transport processes in heterogeneous media, as documented by, for example, Hafskjold et al. (1993); Ikeshoji and Hafskjold (1994); Hafskjold and Ikeshoji (1996).
We use our coarse-graining procedure to obtain the effective pressure of the REV. Its negative gradient is the driving force for fluid flow. The effective pressure of a REV, p , is called the integral pressure. It is in general a combination of pressures, surface tensions, and line tensions. This approach originates from nanothermodynamics, as described by Hill (1994). We believe, however, that this procedure is general for other porous media and fluids. We have shown that the method can replace the use of Young or Young-Laplace's law (Rauter et al. 2020). For a recent formulation of nanothermodynamics, see Bedeaux et al. (2020).
It is our long-range aim to obtain a procedure that provides equations of transport in porous media in general. The aim of this work is not to obtain accurate transport coefficients, but rather to develop a thermodynamic theory for transport in porous media. Focus is therefore given to methods and principles. Not only for isothermal transport of one fluid but also coupled transport due to other driving forces, for example, thermal driving forces (Rauter et al. 1 3 2021, 2021). We start with a single-phase flow, to document the use of new concepts on porous media pressure in a simple way. We have reported equilibrium studies with this coarsegraining procedure previously (Galteland et al. , 2022Erdős et al. 2020;Rauter et al. 2020). We expand these results to non-equilibrium conditions in this work.
This work aims to analyze this simple transport problem and compare its results to expressions that are common in the literature, most importantly Darcy's law in the presence of the Klinkenberg effect Helmig (1997), but also the Kozeny-Carman equation Kozeny (1927);Carman (1937Carman ( , 1956; Berg (2014). Figure 2 illustrates the chosen system, a fcc lattice of solid particles (blue), all with radius R, and a lattice constant a. The porosity is varied by varying a. The lattice is periodic in all directions. The REV is the size of the lattice unit cell ), but here we chose to integrate to a slab of volume V l . The volume V l contains four unit cells and is practical for the purpose. A pressure difference arises using a reflecting particle boundary, which forces a mass flux J m through the cross-sectional area. This boundary-driven method gives minimal disturbance of fluid particles away from the boundary.
In Sect. 2, we recapitulate the necessary thermodynamic equilibrium relations of the pertinent REV, followed by a description of transport, and a procedure to find pressure profiles away from equilibrium. The equilibrium and non-equilibrium simulation procedures that are applied, are described in Sect. 3. The results are discussed in Sect. 4. We give the results needed to determine the local driving force inside the porous medium. The equation of state of the porous medium is central here. The hydraulic conductivity and the permeability are however strong functions of porosity and average fluid mass density. We shall also see that the Klinkenberg correction, as described by Helmig (1997), applies to the fluid permeability. = (237 ± 1)kg∕m 3 . The blue particles represent the solid particles of a fcc lattice, while the red particles represent the fluid particles. There is an integral pressure difference Δp , which drives a mass flux J m from left to right. The slab of volume V l was used as REV. The simulation box had periodic boundaries in all directions (black dashed line), apart from particles crossing the x-boundary from left to right. Here the reflecting particle method was applied, see Sect. 3 for details of this boundary condition (Li et al. 1998).

The Grand Potential and the Pressure
We have chosen to describe the state of a porous medium using a basis set of thermodynamic variables that are additive in the sense that each coarse-grained variable in the REV is a sum of contributions from all bulk phases, surfaces, and possible contact lines (Kjelstrup et al. 2018Bedeaux and Kjelstrup 2021). The REV is an open thermodynamic system, such that mass and heat can be transported in and out of a REV from the surrounding REVs. In the equilibrium situation, the REV is in equilibrium with a bulk fluid with the same temperature and fluid chemical potential. The control variables of the REV are temperature, volume, and fluid chemical potential. The REV of interest belongs therefore to a grand canonical ensemble. The REV variables are controlled by the environment ). The REV is in thermodynamic equilibrium when it has the same temperature and fluid chemical potential as its surroundings. In non-equilibrium conditions, we assume local equilibrium of the REV. In local equilibrium, the temperature and fluid chemical potential are given locally but allowed to vary along the direction of transport and as a function of time. We recapitulate the thermodynamic description of a REV of this kind at equilibrium before we define the situation in the presence of flow.
In nanothermodynamic theory, the grand potential Υ of the REV is given by minus the so-called integral pressure of the REV times the volume. The adjective "integral" was coined by Hill (1994) to reflect the fact that it is an integrated property. Conversely, the commonly used pressure was called the differential pressure. We can regard the integral pressure as the effective pressure of the REV. The grand potential is In this expression, p and V are the integral pressure and the volume of the REV, respectively. It can be understood as the defining equation for p . The grand potential has its basis in statistical mechanics where Ξ is the grand canonical partition function, k B is Boltzmann's constant, and T is the temperature. This basis explains why Υ has additive contributions. Weakly coupled sub-systems will add to Υ , and result in a product of partition functions. This expression explains also why the REV needs to include all possible micro-states of the system.
A more colloquial name of the grand potential is the compressional energy, since it is a product of pressure and volume, which is related to work. In the present case, the grand potential has several additive contributions. We obtain The contributions are from the fluid and solid phases, and the fluid-solid surface. The symbols p f , p s are the integral pressure of the fluid and solid, respectively, and ̂ is the integral fluid-solid surface tension. The surface energies are more significant because of the fluid confinement. For large fluid volumes, the surface energy may be neglected. However, when Young-Laplace's and Young's equations are significant for multi-phase flow the surface energies in this equation are too. The volume of the fluid and solid is V f and V s , respectively, and A is the surface area between the solid and the fluid. The sum of the fluid and (1) Υ = −pV. (2) solid volumes is equal to the REV volume V = V f + V s , and the porosity is = V f ∕V . It follows that the fluid mass density of the REV is where M f is the mass of the fluid in the REV. This fluid mass density follows the coarsegraining procedure. This density is not the internal density given by M f ∕V f . The individual contributions to the integral pressure are not needed, and we only need the total sum. The integral pressure is obtained by calculating an equation of state in equilibrium conditions, which is applied to non-equilibrium conditions. The equation of state for the porous medium must as all equations of state be found from experiments or simulations. In the present work, the last method gives the equation of state for the present medium The functional dependence of p on the variables ( f , ) will here be investigated by simulations at isothermal conditions. The value of p shall also be found by dividing both sides of Eq. 3 by the REV volume. This gives In earlier work we determined p s from known values of A, V and ).
Here we shall use Eq. 6 to compute p and compare the result to the result from Eq. 5. Both procedures are completely general. In Eq. 6 we apply an independent computation of the single parameters. The two expressions will be shown here to give the same result.
The general relation between the differential and integral pressures is given by This shows that p or p can enter the equation of state for the porous medium, see Eq. 5. Because of the lattice symmetry, one unit cell is a large enough REV such that thermodynamic properties do not depend on the size of the REV. The integral pressure therefore will not depend on the REV volume, giving the special condition p =p for the REV. Any block of adjacent unit cells will therefore also give the same value of p , see Fig. 1. It follows that the integral pressure in this lattice is independent of the size of the REV. This is only the case for this special system, and not in general. To be general, we will keep the integral pressure p. Such an equality does not apply for p f , p s and ̂ , which in general differ from p f , p s and . Only if the spheres are far apart, p f = p f , and if the spheres are large enough, also ̂= . This was the special case considered before . In that work, we required that the fluid volume was so large that the fluid pressure was equal to the bulk pressure and that the fluid-solid curvature was so small that the surface tension did not depend on it. By using such requirements, we cannot take into account any disjoining pressure (Israelachvili 2015) or Tolman length (Tolman 1949). The procedure in this work is more general and can describe such effects and other capillary effects.
The presence of additional fluid phases would contribute by similar terms to Eq. 6. The saturation will appear as a variable, and there are contributions from three-phase contact lines. Equation 6 is an alternative definition of the integral pressure of the REV. When there is equilibrium at the boundary between the porous medium and the bulk phase surrounding the porous medium, we have the condition This condition was recently used to determine the solid integral pressure p s (Galteland et al. 2022). We shall here use the relation 8 to determine the REV integral pressure in the equation of state (Eq. 5). The equation of state is next applied to non-equilibrium conditions assuming local equilibrium. This procedure is possible for any geometry.

The Effective Pressure Gradient
The Gibbs equation can be written in terms of coarse-grained variables of the REV, as described in the previous section. By introducing the entropy, mass, and energy balance equations, we can then obtain the entropy production of the REV. In the present isothermal single-phase case, the entropy production has only one flux-force product. The flux conjugate to the negative pressure gradient is the volume flux, J V This expression gives a thermodynamic basis for Darcy's law, which is a locally linear relationship between the volume flux and the driving force valid for the porous medium. According to non-equilibrium thermodynamics, this equation does not assume laminar flow conditions. The only assumption is that there is local equilibrium and that the fluxes are linear combinations of the forces. These assumptions have been shown to hold for coupled heat and mass transport through liquid-vapor interfaces (Røsjorde et al. 2000) and membranes (Inzoli et al. 2009). The conductivity coefficient (the hydraulic conductivity), l, is a function of state variables. The driving force is not necessarily constant along the medium. We need the local value of p∕ x and J V to find l. When the porous medium is in contact with bulk fluids at both ends, the boundary condition will always give Δp = Δp b . In the experimental situation, the gradient is the pressure difference divided by the length L of the medium and the conductivity coefficient is the coefficient l obtained from The permeability k enters Darcy's law on the form where b is the shear viscosity of the bulk fluid, and The volume flux is related to the mass flux by the fluid mass density, We shall determine the coefficients l and k from the last two equations, with information of the mass flux J m , fluid mass density of the REV, f , and integral pressure gradient p∕ x . Their dependence on the fluid mass density f and the porosity is of interest. The original form of Darcy's law was obtained experimentally for the laminar flow of a single-phase fluid through a porous medium. It has also been given a mechanical basis (Helmig 1997). The Reynolds number, used to characterize flow patterns, is here calculated from Re = 2RJ m ∕ b , where J m is the mass flux, and 2R is the diameter of the solid particles. For Reynolds numbers considered in this work, the flow is laminar (Helmig 1997).
A question of general interest is how the permeability of a porous medium, as used in Darcy's law, can be related to the properties of the porous medium. A well-known approach for a bed of granular solids is given by the Kozeny-Carman equation. As the name suggests, it was first derived by Kozeny (1927), assuming that the fluid in the porous medium could be described as contained in a set of non-interfering parallel channels with the same internal surface and pore volume, as the medium itself. The solid phase was a packed bed. The general equation was written as Kozeny (1927); Carman (1937Carman ( , 1956; Berg (2014) where is the porosity, c 0 denotes the Kozeny constant, is the tortuosity, and A/V is the surface-to-volume ratio. The tortuosity is the average length a particle travels across the medium divided by the medium length L. For monodisperse spheres of diameter d the permeability simplifies to Klinkenberg gave a correction for low-density compressible fluids with basis in fluid slippage at the wall (Klinkenberg 1941). When this correction is included, the effective permeability, k, is written as where k 0 is the absolute permeability and b is the Klinkenberg constant (Helmig 1997). The correction was found to be relevant for permeabilities k < 10 −13 m 2 (Baehr and Hult 1991).

The Porous Medium
Systems of a single-phase and single-component fluid were simulated in a porous structure using molecular dynamics simulations with LAMMPS Thompson et al. (2021). The porous solid structure was composed of solid spherical particles in a face-centered (16) k = k 0 1 + b p cubic (fcc) lattice. The fluid particles were free to move, while the solid particles were immovable. As a consequence, the porous medium was non-deformable. The mass of the fluid particles was m, while the mass of the solid particle was not defined since they were immovable. The radius of the solid particles was R and the lattice constant was a.
The particles interacted with the Lennard-Jones/spline potential, which is a pairwise interaction potential that models non-bonded neutral atom interactions. It is, for example, an accurate model for noble gasses and methane (Steele 1974;Miyahara and Gubbins 1997;McGaughey and Kaviany 2004;Rutkai et al. 2017). See Hafskjold et al. (2019) and Kristiansen (2020) for details on the thermodynamic and transport properties for this potential. The potential is where r = ‖r j − r i ‖ was the distance between particle i and j. The interaction strength was characterized by the depth of the interaction potential , and the particle diameter was characterized by the distance . The hard-core diameter of the particles was d. The parameters a, b, r s and r c were set such that the potential and its derivative were continuous at r s and r c .
The microscopic mechanical pressure tensor is ambiguous. The virial pressure due to pair interactions can be distributed along any line connecting particle pairs (Harasima 1958;Schofield and Henderson 1982;Long et al. 2011;van Dijk 2020;Long et al. 2020). The mechanical pressure is well-defined in the bulk fluid. Because of this, we have used the equation of state approach. The mechanical pressure is calculated in the bulk fluid which is in equilibrium with a porous medium. We assume the integral pressure of the porous medium to be equal to the mean of the diagonal components of the mechanical pressure of the bulk fluid.
The simulation data were related to SI units by using parameters , , and the fluid particle mass m for methane Steele (1974); Miyahara and Gubbins (1997). The aim of this is not to give accurate data for methane transport, but rather to present more identifiable values. The parameters were in terms of SI units The parameters ff = fs = and ff = fs = were equal for the fluid-fluid and solidfluid interactions, while the hard-core diameter was zero for the fluid-fluid interactions d ff = 0 and d fs = 4.5 ≈ 1.7 nm for the fluid-solid interactions. The solid did not interact with other solid particles. The radius of the solid was constant and defined to be R ∶= d fs + fs ∕2 = 5 ≈ 1.9 nm , at which point the potential energy of the fluid-solid interaction is zero. The fluid-solid slip conditions were not precisely defined, but depended on the shape of the fluid-solid surface and fluid-solid interactions. The fluid-solid surface was completely smooth.
The lattice constant was varied from a = 2 √ 2R ≈ 5.4 nm to a = 40 ≈ 15.2 nm , where the lower limit is where the surface of the solid particles are in contact. The porosity and surface-to-volume ratio of the porous structure were (18) ∕k B = 148.1 K, = 0.381 nm, and m = 2.661 × 10 −26 kg.
Consequently, the porosity varied from 0.26 to 0.97 and the surface-to-volume ratio varied between approximately 0.05nm −1 and 1.17nm −1 . The fcc lattice has four octahedral and eight tetrahedral voids. The octahedral voids are the largest and can accommodate a sphere of radius a∕2 − R , which varies from 0.8 nm to 5.7 nm in the cases studied here.

Shear Viscosity
The shear viscosity of the bulk fluid was calculated with the OCTP plugin (Jamali et al. 2019) from the time-integral of the auto-correlation function of the off-diagonal components of the Cartesian mechanical pressure tensor, The system consisted of 32000 fluid particles and was initialized at temperature T = 296.2 K and bulk mass density in the range 4.8 to 385kg∕m 3 . The system was equilibrated with an NVT-ensemble for 2 × 10 6 steps. After that, an additional 2 × 10 6 steps were run in the NVE ensemble to collect data and compute the shear viscosity. The time step was t = 1.375 fs . The statistics were improved by running 30 independent simulations. The shear viscosity is shown as a function of the fluid mass density in Fig. 3.

Equilibrium Conditions
The systems were initialized by creating a fcc unit cell of solid particles with varying lattice constant a in a cubic simulation box of side lengths a. In addition, a bulk fluid without solid particles was simulated. The boundaries of the simulation box were periodic. The temperature of the fluid was controlled with a Nosé-Hoover-type thermostat to be constant and equal to T = 2 ∕k B = 296.2 K (Shinoda et al. 2004), which is in the Fluid particles were inserted and removed from the porous structure using the grand canonical Monte Carlo (Frenkel and Smit 2002). This was done to generate initial configurations with varying lattice constants that were in equilibrium with each other. Fluid particles were inserted into the pores with an acceptance probability and a random fluid particle was removed from the simulation box with an acceptance probability where = 1∕k B T is the thermodynamic beta, Λ = √ h 2 ∕(2 mk B T) is the de Broglie thermal wavelength, h is the Planck constant, N f the number of fluid particles, and ΔE p is the potential energy difference of the system, when inserting or removing a fluid particle. The fluid chemical potential varied in the interval f ∈ [−10, 10] , which resulted in a bulk fluid mass density variation from b = (3.3 ± 0.2) kg∕m 3 to (409.9 ± 0.5) kg∕m 3 . The simulations were run until the average number of fluid particles reached a constant value. For the largest porosities and chemical potentials, this took up to 9 × 10 6 steps. The time step was t = 2.75 fs.  Each illustration is scaled such that the figure sizes become the same. In the simulations, the radius of the blue solid particles was constant in all simulations. The simulated system was visualized with OVITO (Stukowski 2009) scaled such that the figure sizes become the same. In the simulations, the radius of all blue solid particles was the same.

Non-equilibrium Conditions
The unit cells were replicated twice in the y-and z-directions and ten times in the x-direction in these studies. The final side lengths of the simulation box in non-equilibrium conditions were (10, 2, 2)a, where a is the lattice constant. The grand canonical Monte Carlo fluid particle insertions and removals from the porous medium were stopped, while the temperature was controlled with a Nosé-Hoover-type thermostat (Shinoda et al. 2004). However, the temperature was also controlled separately in each REV (volume V l ) to ensure that there was no temperature gradient. The reflecting particle method (RPM) (Li et al. 1998) was applied to the x-boundary of the simulation box to induce an integral pressure gradient. The RPM allowed particles to cross the periodic x-boundary from right to left with a probability 1 − q and be reflected with a probability q ∈ [0, 1] . A value q = 0 entailed that no fluid particles were reflected, while a value q = 1 entailed that all fluid particles were reflected when attempting to cross the x-boundary from right to left. The fluid particles were never reflected when crossing the x-boundary from left to right. The pressure gradient was controlled by varying the probability q. This is a boundary-driven non-equilibrium molecular dynamics method to induce a pressure gradient that has minimal disturbance on the fluid particles away from the boundary. The simulations were run until they had reached steady state, where the mass flux was constant along the x-axis. For the longest simulations, this took up to 8 × 10 6 steps. The time step was t = 2.75fs.
The mass fluxes and integral pressures were calculated in layers l of volume V l along the x-axis. The volume of each layer was equal to the lattice constant a and spanned the simulation box in the y-and z-directions. The side lengths of the layers were (a, 2a, 2a) in the x-, y-and z-directions, respectively. The volume of each layer (the REV) was consequently V l = 4a 3 . Each layer contains four unit cells. The mass flux through layer l was calculated as where m i is the mass of fluid particle i, v i,x is the velocity of fluid particle i in the x-direction, and ⟨v i,x ⟩ is the average particle velocity in the x-direction. The sum is over all particles in the layer V l . The integral pressure of each layer was calculated by relating the fluid mass density f ,l to the integral pressure calculated in equilibrium.
With this procedure, we were able to determine the transport coefficients l and k. The procedure has been used earlier by us to establish the same gradient in integral pressure (Galteland et al. 2022); however, the permeability calculations are done for the first time here.

Results and Discussion
The results of this work are shown in Figs. 5,6,7,8,9,10,11,12,13 and 14. The results are presented and discussed with reference to the simulation procedure (Sect. 3.3) and to the Theory (Sect. 2).

Equilibrium Conditions: The Equation of State
The fluid mass density in the REV is shown as a function of fluid chemical potential in Fig. 5. The chemical potential was set by the environment. The fluid mass density is shown for a varying porosity. We see that the fluid mass density in the REV increases The fluid mass density of the REV differs from the mass density outside the porous medium. But, the fluid mass density of the REV, corrected by the porosity will also . This indicates that thermodynamic properties change upon fluid confinement. In the low porosity case, the effect is large. For a fluid chemical potential near 0 J, the density f changes by a factor of 50, as the porosity changes from the densest packing of spheres to bulk fluid. Every single curve deviates from the bulk vapor mass density when the fluid chemical potential is The dependence of the integral pressure on the fluid mass density for the various porosities is presented in Fig. 6. The rise in the integral pressure from a few bars to (3184 ± 27)bar is shown for all porosities. The slope of the curves increases with decreasing porosity. The bulk isotherm is again indicated by black crosses and gives the lower limit of the variation. Again, the curves approach the bulk value as the porosity increases, as expected. The curves apply for a temperature T = 2.0 ∕k B = 296.2 K . The isothermal compressibility T = −1 f f ∕p T of the bulk fluid varied from (0.20 ± 0.03)bar −1 to (0.753 ± 0.001) × 10 −7 bar −1 . See the supplementary information for details.
The mean free path of the bulk vapor densities can be estimated with the kinetic theory of gases with the equation = ( √ 2 2 n b ) −1 , where n b is the number density of fluid particles in the bulk phase. For the lowest bulk mass density, b = (3.3 ± 0.2)kg m −1 the mean free path is estimated to be = (12.7 ± 0.9) nm . The octahedral voids in the fcc lattice can accommodate spheres with diameter 1.6 nm to 11.4 nm for the varying porosities, in which case the system is in the Knudsen flow regime. For a bulk mass density b = (24.2 ± 0.6)kg m −1 , the mean free path is estimated to be = (1.71 ± 0.05)nm . The system will be in the Knudsen flow regime for the lowest densities.
The curves obtained for the equilibrium condition Eq. 6 can be regarded as an equation of state relating the temperature, porosity, and fluid mass density of the REV to the integral pressure, see equation Eq. 5. we shall use the sets of relations in the same manner as an equation of state. Once we know the fluid mass density of the REV, we know also the chemical potential that controls it from Fig. 5, and therefore also the corresponding integral pressure from Fig. 6.

Non-equilibrium Conditions
By applying the reflective boundary method we generated a mass flux, a fluid mass density difference, and a difference in integral pressure difference across the porous medium. By increasing the reflecting probability q we increased the gradients. In this manner, we varied the integral pressure gradient between approximately Δp∕L = −7 bar m −1 and 20000 bar m −1 . The corresponding Reynolds numbers varied from approximately zero up to Re = 3.54 ± 0.08 . Gradients that are generated in non-equilibrium molecular The fluid mass density profiles across the porous medium are shown in Fig. 7 for a porosity ≈ 0.87 , and an average fluid mass density of the REV f = (351.3 ± 0.7)kg∕m 3 . The Reynolds numbers varied from Re = 0.02 ± 0.01 to Re = 0.67 ± 0.02 . From the fluid mass density profiles, we obtained the integral pressure profiles. Examples are shown in Fig. 8. The integral pressure gradients were calculated from the integral pressure profiles. The gradient in integral pressure was in good approximation constant. This was not expected and is also not needed in the data reduction procedures. Because the mass flux is constant at a steady state (mass conservation), and the driving force is approximately constant, it follows that f l and f k∕ were constant across the porous medium. The conductivity and the permeability are not necessarily constant in a porous medium of varying porosity. There is not yet much experience with the integral pressure of porous media reported in the literature. It is therefore appealing to examine the contributions to the integral pressure from the bulk phases and surface, and test our way to compute p from the equation of state, Eq. 5. By assuming p f = p f , p = p f and ̂= , the individual contributions can be calculated. The outcome is illustrated in Fig. 9, again for a porosity ≈ 0.87 , average fluid mass density b = (351.3 ± 0.7)kg∕m 3 . The Reynolds number was now Re = 0.67 ± 0.02 . The assumptions hold for large porosities when the disjoining pressure is zero, and when the surface tension is independent of the fluid-solid surface curvature. These assumptions are not necessary for the continued analysis of this work, they are just chosen to illustrate a case with individual contributions. The surface tension was calculated from the spherical mechanical pressure tensor. The value increases monotonically with increasing density between (0.079 ± 0.009)mN m −1 and (34 ± 1)mN m −1 . See the supplementary information for more details on the surface tension.
The integral solid pressure was calculated from the integral pressure and surface tension where A is the fluid-solid surface area and V s is the volume of the solid phase, see Eq. 3. In Fig. 9, we see the local contributions from the gradient in fluid pressure ∇p f = (−8629 ± 106)bar m −1 , in integral solid pressure (1 − )∇p s = (−9965 ± 215)bar m −1 , and in surface tension (A∕V)Δ = (370 ± 214)bar m −1 . The single parts sum to the total value shown in the figure. This value of p was also determined from the equation of state, Eq. 5. Within the accuracy of the calculation, we confirmed our hypothesis that the two routes give the same result.
A control of isothermal conditions was carried out. We know that a pressure difference may generate a temperature difference, or vice versa ). Constant temperature is, therefore, a condition for the single flux-force product of the entropy production, and consequently for Darcy's law's applicability. We, therefore, confirmed that the temperature was constant across the system. Figure 10 shows plots that are used to determine the conductivity and permeability in Darcy's law. Rather than using the volume flux, we have used the mass flux on the ordinate axis, as this flux, but not the volume flux is constant across the porous media in a steady state. The figure illustrates varying porosities ≈ [0.26, 0.61] and the average bulk fluid mass density was b = (108 ± 1)kg∕m 3 in these plots.
(24) p s = p f + A∕V s We see that the mass fluxes in all cases can be regarded as linear functions of the negative integral pressure gradient; the dashed lines are linear fits to the calculated mass fluxes. Also observed is that the fluxes intersect the origin of the axes. This behavior is compatible with a linear theory like NET. There is no threshold for transport at low-pressure gradients. Figure 8 gives the integral pressure gradient inside the porous medium. The integral pressure gradient from this figure together with the constant mass flux was used to determine the conductivity of the porous medium. The conductivity is plotted in Fig. 11 as a function of the factor 3 ∕(1 − ) 2 in the Kozeny-Carman equation. There was no convincing relationship between the variables.

Conductivity and Permeability
It is more interesting to consider the permeability k. Its value was computed from the conductivity l and the shear viscosity b , and the results are shown as a function of the inverse average integral pressure (which here happens to be the same as the inverse differential pressure) in Fig. 12. We see that the permeability increases with a constant slope for large values of the inverse average integral pressure (low densities), where the fluid is highly compressible. For small values (high densities), the slope decreases as the fluid become less compressible. This behavior is expected from the Klinkenberg correction (Klinkenberg 1941;Baehr and Hult 1991). The results at the low pressure were therefore fitted to a straight line and extrapolated to the high-pressure end to obtain the absolute permeability, k 0 .
The dependence of the permeability on other lattice parameters is shown in the log-plot of the permeability as a function of the porosity in Fig. 13. A systematic variation is demonstrated over a change in four orders of magnitude of k. The absolute permeability was calculated using the Klinkenberg correction formula (in the range where this correction applies), by plotting the permeability vs 1/p and extrapolating to very large pressures. The absolute permeability is shown as a black dashed line. We see that the absolute permeabilities form a lower limit for the family of calculated permeabilities. On the other hand, the Kozeny-Carman equation with tortuosity = 1 , shown as a grey dash-dotted line, does not bring out any new physical insight. Figure 14 shows the Klinkenberg coefficient b as a function of porosity . The coefficient is porosity dependent and significant for the whole range of porosities. To a good approximation, the dependence is linear for porosities below 0.6. The coefficient can be understood from the lack of interaction between solid spheres and fluid particles on the particle level, which here is consistent with slippage.

Conclusion and Perspectives
In this work, we have developed a method of investigating transport process with a basis in nanothermodynamics and non-equilibrium thermodynamics and shown how this can be used to obtain transport properties. We have obtained the integral pressure of the REV in non-equilibrium conditions by using an equation of state constructed in equilibrium conditions. To construct the equation of state we have assumed the integral pressure to be equal in the REV and a bulk fluid with the same temperature and fluid chemical potential. The REV was constructed from additive thermodynamic variables, such that the Gibbs equation and in turn the entropy production and flux-force equations could be derived. The method was applied to single-phase fluid flow in an isothermal medium of varying porosity and fluid mass density under laminar flow conditions. The hydraulic conductivity and the permeability were shown to vary with porosity and fluid mass density, and give values that are typical in the literature. The systems studied in this work were in steady-state. This does not pose any limitation on the method. The procedure can be used to describe transients. It is interesting that the system supported the behavior behind the Klinkenberg effect, and that we obtained variables of a size compatible with this effect, e.g. k < 10 −13 m 2 .
The procedure has its thermodynamic basis in standard non-equilibrium thermodynamics as combined with Hill's method of nanothermodynamics to describe the confined fluid. This method should now be tested with a two-phase flow, to help solve problems stated in the literature on the upscaling problem; i.e. on how we can properly describe the porous medium microstates, that are the origin of Darcy scale behavior. For instance, we can define the effective pressure without consideration of the capillary pressure. For twophase flow, the REV will be larger than a unit cell to include the statistical variation of the two-phases.
A particular symmetric lattice was chosen to illustrate the derivations. This should also not be regarded as a limitation. The method could be accommodated to deal with pore distributions. It can be applied to highly confined systems, where the disjoining pressure is significant. The central point is the proper construction of the grand potential.