# Local thermodynamic description of isothermal single-phase flow in porous
media

Olav Galteland  PoreLab, Department of Chemistry, Norwegian University of
Science and Technology  Corresponding author<EMAIL_ADDRESS>Michael
T. Rauter  PoreLab, Department of Chemistry, Norwegian University of Science
and Technology  Mina S. Bratvold  PoreLab, Department of Chemistry, Norwegian
University of Science and Technology  Thuat T. Trinh  PoreLab, Department of
Chemistry, Norwegian University of Science and Technology  Dick Bedeaux
PoreLab, Department of Chemistry, Norwegian University of Science and
Technology  Signe Kjelstrup  PoreLab, Department of Chemistry, Norwegian
University of Science and Technology

###### Abstract

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\text{\times}{10}^{-20}\text{\,}\mathrm{m}^{2}$ and
$5.5\text{\times}{10}^{-17}\text{\,}\mathrm{m}^{2}$, 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.

Keywords: nanothermodynamics, Hill’s thermodynamics of small systems, non-
equilibrium thermodynamics, non-equilibrium molecular dynamics, nanoporous
media, compressible flow, representative elementary volume

##  1 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 [1]. We also need to describe the selective transport across the porous
separators in batteries and fuel cells [2]. We are aiming for a description
that reflects the underlying properties of the single pores.

Much of the fundamental work on flow in porous media has been done for the
pore-scale, see Helmig for a pedagogical presentation [3]. Recent developments
in imaging techniques and computer capabilities have improved our
understanding of the physics at pore-scale considerably [4]. 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 the fluid. We shall use additive variables to define the
REV, similar to Whitaker [5]. This definition of the REV differs from the one
suggested by Nordahl and Ringrose [6] who used a constant permeability as
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, $\hat{p}(T,\rho_{f},\phi)$. This can be regarded as finding the
equation of state of the REV. Up to this point, we have used REV densities
much like Whitaker [5].

Figure 1: The representative elementary volume (REV) for a single-phase fluid
in a face-centered 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 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 [7]. 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 [8,
9, 10].

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 [11, 12, 10]. We have
derived the Gibbs equation for the REV, and in the analysis of transport, we
have assumed validity of Gibbs equation. This gives the entropy production for
the REV and a thermodynamic basis for Darcy’s law in 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 which
is standard in non-equilibrium thermodynamics.

An excellent tool to analyze transport in porous media is non-equilibrium
molecular dynamics (NEMD) [13]. 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 becomes limited as NEMD is a computationally
expensive technique. Here, we will use NEMD to simulate the flow of methane-
like molecules in a face-centered cubic (fcc) lattice made up of spherical
solid particles, 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, Ikeshoji and Hafskjold [14, 15, 16].

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, $\hat{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 [17]. We
believe, however, that this procedure for other porous media and other fluids.
We have shown that the method can replace the use of Young or Young-Laplace’s
law [9]. For a recent formulation of nanothermodynamics, see Bedeaux,
Kjelstrup, and Schnell [18].

It is our long-range aim to obtain a procedure that provides equations of
transport in porous media in general. Not only for isothermal transport of one
fluid but also coupled transport due to other driving forces, for example,
thermal driving forces [19, 20]. 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 coarse-graining procedure
previously [8, 21, 9, 22, 23, 24]. We expand these results to non-equilibrium
conditions in this work.

Figure 2: A visualization of a simulation, with porosity $\phi\approx 0.87$
and average bulk fluid mass density $\rho_{b}=$237\pm
1.\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{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
$\Delta\hat{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 section 3 for details of this boundary condition [25]. The
simulated system was visualized with OVITO [26].

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 [3], but also the
Kozeny-Carman equation [27, 28, 29, 30].

Fig. 2 illustrates the chosen system, a fcc lattice of solid particles (blue),
all with radius $R$, and a lattice constant $a$. The porosity $\phi$ is varied
by varying $a$. The lattice is periodic in all directions. The REV is the size
of the lattice unit cell [8], 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 Section 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
Section 3. The results are discussed in Section 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 [3], applies to the fluid permeability.

##  2 The thermodynamic variables of the REV

###  2.1 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 [11, 12, 10]. The REV variables are
controlled by the environment [18]. In the present case, the control variables
are temperature, volume, and fluid chemical potential. The REV of interest
belongs therefore to a grand canonical ensemble. 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 $\Upsilon$ 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 [17] to reflect the fact that it
is an integrated property. Conversely, the normal pressure was called the
differential pressure. We can regard the integral pressure as the effective
pressure of the REV. The grand potential is,

$\Upsilon=-\hat{p}V.$ (1)

In this expression, $\hat{p}$ and $V$ are the integral pressure and the volume
of the REV, respectively. It can be understood as the defining equation for
$\hat{p}$. The grand potential has its basis in statistical mechanics

$\Upsilon=-k_{\text{B}}T\ln\Xi,$ (2)

where $\Xi$ is the grand canonical partition function, $k_{\text{B}}$ is
Boltzmann’s constant, and $T$ is the temperature. This basis explains why
$\Upsilon$ has additive contributions. Weakly coupled sub-systems will add to
$\Upsilon$, 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

$\hat{p}V=\hat{p}_{f}V_{f}+\hat{p}_{s}V_{s}-\hat{\gamma}A$ (3)

The contributions are from the fluid and solid phases, and the fluid-solid
surface. The symbols $\hat{p}_{f}$, $\hat{p}_{s}$ are the integral pressure of
the fluid and solid, respectively, and $\hat{\gamma}$ 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 are $V_{f}$ and $V_{s}$, respectively, and $A$ is the
surface area between the solid and the fluid. The sum of the fluid and solid
volumes is equal to the REV volume $V=V_{f}+V_{s}$, and the porosity is
$\phi=V_{f}/V$. It follows that the fluid mass density of the REV is

$\rho_{f}=\frac{M_{f}}{V_{f}+V_{s}}$ (4)

where $M_{f}$ is the mass of the fluid in the REV. This fluid mass density
follows the coarse-graining procedure. This density is not the internal
density given by $M_{f}/V_{f}$. The individual contributions to the integral
pressure are not needed, 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

$\hat{p}=\hat{p}(T,\rho_{f},\phi).$ (5)

The functional dependence of $\hat{p}$ on the variables ($\rho_{f},\phi$) will
here be investigated by simulations at isothermal conditions. The value of
$\hat{p}$ shall also be found by dividing both sides of Eq. 3 by the REV
volume. This gives

$\hat{p}=\hat{p}_{f}\phi+\hat{p}_{s}(1-\phi)-\hat{\gamma}A/V.$ (6)

In earlier work we determined $\hat{p}_{s}$ from known values of $A,V$ and
$\phi$ [8]. Here we shall use Eq. 6 to compute $\hat{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

$p=\left[\frac{\partial(\hat{p}V)}{\partial
V}\right]_{T,\mu}=\hat{p}+V\left(\frac{\partial\hat{p}}{\partial
V}\right)_{T,\mu}$ (7)

This shows that $p$ or $\hat{p}$ can enter the equation of state for the
porous medium, see Eq. 5. The fcc lattice used in the present model poses a
special condition on the two pressures. Because of the lattice symmetry, the
integral pressure will not depend on the REV volume, giving the special
condition $p=\hat{p}$ for the REV. Any block of adjacent unit cells will
therefore also give the same value of $\hat{p}$, see Fig. 1. It follows that
the integral pressure in this lattice is independent of the size of the REV.
As a consequence, the differential pressure of the REV is identical to the
integral pressure, $p=\hat{p}$ for the REV. This is only the case for this
special system, and not in general. To be general, we will keep the integral
pressure $\hat{p}$.

Such an equality does not apply for $\hat{p}_{f}$, $\hat{p}_{s}$ and
$\hat{\gamma}$, which in general differ from $p_{f},p_{s}$ and $\gamma$. Only
if the spheres are far apart, $\hat{p}_{f}=p_{f}$, and if the spheres are
large enough, also $\hat{\gamma}=\gamma$. This was the special case considered
before [8]. 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 [31] or
Tolman length [32]. 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. Eq. 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

$\hat{p}=p_{b}.$ (8)

This condition was recently used to determine the solid integral pressure
$\hat{p}_{s}$ [22]. 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.

###  2.2 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}$

$J_{V}=-l\frac{\partial\hat{p}}{\partial x}.$ (9)

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 [33] and membranes [34]. 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 $\partial\hat{p}/\partial 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 $\Delta\hat{p}=\Delta 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

$J_{V}=-l\frac{\Delta\hat{p}}{L}=-l\frac{\Delta p_{b}}{L}.$ (10)

The permeability $k$ enters Darcy’s law on the form

$J_{V}=-\frac{k}{\eta_{b}}\frac{\Delta\hat{p}}{L}=-\frac{k}{\eta_{b}}\frac{\Delta{p}_{b}}{L},$
(11)

where $\eta_{b}$ is the shear viscosity of the bulk fluid, and

$l=\frac{k}{\eta_{b}}$ (12)

The volume flux is related to the mass flux by the fluid mass density,
$J_{V}=J_{m}/\rho_{f}$. While $J_{V}$ varies, $J_{m}$ is constant in the
steady state,

$J_{m}={J_{V}}{\rho_{f}}=-\rho_{f}\frac{k}{\eta_{b}}\frac{\partial\hat{p}}{\partial
x}.$ (13)

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,
$\rho_{f}$, and integral pressure gradient $\partial\hat{p}/\partial x$. Their
dependence on the fluid mass density $\rho_{f}$ and the porosity $\phi$ 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 [3]. The Reynolds number, used to characterize flow patterns,
is here calculated from $\text{Re}=2RJ_{m}/\eta_{b}$, where $J_{m}$ is the
mass flux, and $2R$ the diameter of the solid particles. For Reynolds numbers
considered in this work, the flow is laminar [3].

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 [27], 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 [27, 28, 29, 30]

$J_{V}=-\frac{\tau^{2}V^{2}\phi^{3}}{c_{0}A^{2}}\frac{1}{\eta_{b}}\frac{\Delta
p}{L}$ (14)

where $\phi$ is the porosity, $c_{0}$ denotes the Kozeny constant, $\tau$ 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

$k=\frac{\tau^{2}d^{2}}{c_{0}}\frac{\phi^{3}}{(1-\phi)^{2}}.$ (15)

Klinkenberg gave a correction for low-density compressible fluids with basis
in fluid slippage at the wall [35]. When this correction is included, the
effective permeability, $k$, is written as

$k=k_{0}\left(1+\frac{b}{p}\right)$ (16)

where $k_{0}$ is the absolute permeability and $b$ is the Klinkenberg constant
[3]. The correction was found to be relevant for permeabilities $k$¡ 10 -13 m2
[36].

##  3 Molecular simulations

###  3.1 The porous medium

Systems of a single-phase and single-component fluid were simulated in a
porous structure using molecular dynamics simulations with LAMMPS [37]. The
porous solid structure was composed of solid spherical particles in a face-
centered 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
pair-wise interaction potential that models non-bonded neutral atom
interactions. It is for example an accurate model for noble gasses and methane
[38, 39, 40, 41]. See Hafskjold et al. [42] and Kristiansen [43] for details
on the thermodynamic and transport properties for this potential. The
potential is

$u(r)=\begin{cases}\infty&\text{if }r\leq d,\\\
4\epsilon\left[\left(\frac{\sigma}{r-d}\right)^{12}-\left(\frac{\sigma}{r-d}\right)^{6}\right]&\text{if
}d<r\leq r_{s},\\\ a(r-r_{c})^{2}+b(r-r_{c})^{3}&\text{if }r_{s}<r\leq
r_{c},\\\ 0&\text{if }r>r_{c}.\end{cases}$ (17)

Where $r=\|\bm{r}_{j}-\bm{r}_{i}\|$ was the distance between particle $i$ and
$j$. The interaction strength was characterized by the depth of the
interaction potential $\epsilon$, and the particle diameter was characterized
by the distance $\sigma$. 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 simulation data were related to SI units by using parameters $\epsilon$,
$\sigma$, and the fluid particle mass $m$ for methane [38, 39]. 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

$\epsilon/k_{\text{B}}=$148.1\text{\,}\mathrm{K}$,\quad\sigma=$0.381\text{\,}\mathrm{nm}$,\quad\text{and}\quad
m=$2.661\text{\times}{10}^{-26}\text{\,}\mathrm{kg}$.$ (18)

The parameters $\epsilon_{\text{ff}}=\epsilon_{\text{fs}}=\epsilon$ and
$\sigma_{\text{ff}}=\sigma_{\text{fs}}=\sigma$ were equal for the fluid-fluid
and solid-fluid interactions, while the hard-core diameter was zero for the
fluid-fluid interactions $d_{\text{ff}}=0$ and
$d_{\text{fs}}=4.5\sigma\approx$1.7\text{\,}\mathrm{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\vcentcolon=d_{\text{fs}}+\sigma_{\text{fs}}/2=5\sigma\approx$1.9\text{\,}\mathrm{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\sqrt{2}R\approx$5.4\text{\,}\mathrm{nm}$$ to
$a=40\sigma\approx$15.2\text{\,}\mathrm{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

$\phi=1-\frac{16\pi R^{3}}{3a^{3}}\quad\text{and}\quad A/V=\frac{16\pi
R^{2}}{a^{3}}.$ (19)

Consequently, the porosity varied from $0.26$ to $0.97$ and the surface-to-
volume ratio varied between approximately $0.05\text{\,}{\mathrm{nm}}^{-1}$
and $1.17\text{\,}{\mathrm{nm}}^{-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\text{\,}\mathrm{nm}$ to $5.7\text{\,}\mathrm{nm}$ in the cases studied
here.

###  3.2 Shear viscosity

The shear viscosity of the bulk fluid was calculated with the OCTP plugin [44]
from the time-integral of the auto-correlation function of the off-diagonal
components of the Cartesian mechanical pressure tensor,

$\eta_{\alpha\beta}=\lim_{t\rightarrow\infty}\frac{1}{2t}\frac{V}{k_{B}T}\left\langle\left(\int_{0}^{t}P_{\alpha\beta}\left(t^{\prime}\right)\mathrm{d}t^{\prime}\right)^{2}\right\rangle.$
(20)

The system consisted of 32000 fluid particles and was initialized at
temperature $T=$296.2\text{\,}\mathrm{K}$$ and bulk mass density in the range
$4.8$ to $385\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$. The
system was equilibrated with an $NVT$-ensemble for $2\times 10^{6}$ steps.
After that, an additional $2\times 10^{6}$ steps were run in the $NVE$
ensemble to collect data and compute the shear viscosity. The time step was
$\delta t=$1.375\text{\,}\mathrm{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.

Figure 3: Shear viscosity as a function of fluid mass density.

###  3.3 Simulation procedure

####  3.3.1 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\epsilon/k_{\text{B}}=$296.2\text{\,}\mathrm{K}$$ [45], which is in the
supercritical region
($T_{c}=0.885\epsilon/k_{\text{B}}\approx$131\text{\,}\mathrm{K}$$ [42]). A
supercritical fluid was investigated to avoid phase separation.

Fluid particles were inserted and removed from the porous structure using the
grand canonical Monte Carlo [46]. 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

$\text{acc}(N_{f}\rightarrow
N_{f}+1)=\min\left\\{1,\quad\frac{V}{\Lambda^{3}(N_{f}+1)}\exp[\beta(\mu_{f}-\Delta
E_{p}]\right\\}$ (21)

and a random fluid particle was removed from the simulation box with an
acceptance probability

$\text{acc}(N_{f}\rightarrow
N_{f}-1)=\min\left\\{1,\quad\frac{\Lambda^{3}(N_{f}+1)}{V}\exp[-\beta(\mu_{f}+\Delta
E_{p}]\right\\}.$ (22)

Where $\beta=1/k_{\text{B}}T$ is the thermodynamic beta,
$\Lambda=\sqrt{h^{2}/(2\pi mk_{\text{B}}T)}$ is the de Broglie thermal
wavelength, $h$ is the Planck constant, $N_{f}$ the number of fluid particles,
and $\Delta 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 $\mu_{f}\in[-10,10]\epsilon$, which resulted in a bulk fluid mass
density variation from $\rho_{b}=$3.3\pm
0.2\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$ to $409.9\pm
0.5\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{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\times
10^{6}$ steps. The time step was $\delta t=$2.75\text{\,}\mathrm{fs}$$.

Fig. 4 illustrates nine visualizations of the equilibrium simulations for
three porosities and three chemical potentials. The figure provides a visual
impression of the range of conditions studied, from closest packing of spheres
and porosity $\phi\approx 0.26$ (bottom row) to an open fcc lattice with
porosity $\phi\approx 0.97$ (top row). The three chemical potentials used in
the snapshots correspond to bulk mass densities $\rho_{b}=[(39.9\pm
0.8),(236\pm 1),(350.9\pm
0.8)]$\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$ from left to
right. The figure parts are scaled such that the figure sizes become the same.
In the simulations, the radius of all blue solid particles was the same.

Figure 4: Nine initialized simulations in equilibrium conditions of varying
fluid chemical potential and porosity. The three chemical potentials used in
the snapshots correspond to bulk mass densities $\rho_{b}=[(39.9\pm
0.8),(236.0\pm 1.1),(350.9\pm
0.8)]$\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$ from left to
right. 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 [26].

####  3.3.2 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 [45]. 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) [25] 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\in[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\times 10^{6}$ steps. The time step was $\delta
t=$2.75\text{\,}\mathrm{fs}$$.

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

$J_{l,m}=\frac{1}{V_{l}}\sum_{i\in V_{l}}m_{i}v_{i,x}=\rho_{f,l}\langle
v_{i,x}\rangle$ (23)

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 $\langle v_{i,x}\rangle$ 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 $\rho_{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 [24], however, the permeability calculations are
done for the first time here.

##  4 Results and discussion

The results of this work are shown in Figs. 5 \- 14. The results are presented
and discussed with reference to the simulation procedure (Section 3.3) and to
the Theory (Section 2).

###  4.1 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 monotonically with the chemical potential,
approaching the value of the adjacent bulk phase (shown by black crosses in
the figure). This approach is as expected.

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 differ from the outside value (not shown in the figure).
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\text{\,}\mathrm{J}$, the density $\rho_{f}$ changes by a factor of 50,
as the porosity changes from densest packing of spheres to bulk fluid. Every
single curve deviates from the bulk vapor mass density when the fluid chemical
potential is around $-{10}^{20}\text{\,}\mathrm{J}$, where the adsorption of
fluid into the pores starts to seriously increase. The numbers refer to the
chosen model; methane.

Figure 5: Fluid mass density as a function of fluid chemical potential for
varying porosities. Bulk fluid values are indicated by black crosses at the
top.

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\pm 27.\text{\,}\mathrm{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\epsilon/k_{\text{B}}=$296.2\text{\,}\mathrm{K}$$.

The mean free path of the bulk vapor densities can be estimated with the
kinetic theory of gases with the equation
$\ell=(\sqrt{2}\pi\sigma^{2}n_{b})^{-1}$, where $n_{b}$ is the number density
of fluid particles in the bulk phase. For the lowest bulk mass density,
$\rho_{b}=$3.3\pm 0.2\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-1}$$ the mean
free path is estimated to be $\ell=$12.7\pm 0.9\text{\,}\mathrm{nm}$$. The
octahedral voids in the fcc lattice can accommodate spheres with diameter
$1.6\text{\,}\mathrm{nm}$ to $11.4\text{\,}\mathrm{nm}$ for the varying
porosities, in which case the system is in the Knudsen flow regime. For a bulk
mass density $\rho_{b}=$24.2\pm
0.6\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-1}$$, the mean free path is
estimated to be $\ell=$1.71\pm 0.05\text{\,}\mathrm{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.

Figure 6: Integral pressure as a function of fluid mass density for varying
porosities, or equation of state for the porous medium. Bulk fluid values are
shown by black crosses at the bottom.

###  4.2 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
$\Delta\hat{p}/L=$-7\text{\,}\mathrm{bar}\text{\,}{\mathrm{\SIUnitSymbolMicro
m}}^{-1}$$ and
$20\,000\text{\,}\mathrm{bar}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$.
The corresponding Reynolds numbers varied from approximately zero up to
$\text{Re}=3.54\pm 0.08$. Gradients that are generated in non-equilibrium
molecular dynamics can be very large because the length of the simulation box
is relatively short and the method can apply strong forces.

Figure 7: Fluid mass density as a function of the $x$-coordinate for porosity
$\phi\approx 0.87$, average fluid mass density $\rho_{b}=$351.3\pm
0.7\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$, and varying
Reynolds numbers.

The fluid mass density profiles across the porous medium are shown in Fig. 7
for a porosity $\phi\approx 0.87$, and an average fluid mass density of the
REV $\rho_{f}=$351.3\pm
0.7\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$. The Reynolds
numbers varied from $\text{Re}=0.02\pm 0.01$ to $\text{Re}=0.67\pm 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 $\rho_{f}l$ and $\rho_{f}k/\eta$ were constant across the porous
medium. The conductivity and the permeability are not necessarily constant in
a porous medium of varying porosity.

Figure 8: Integral pressure as a function of the $x$-coordinate for porosity
$\phi\approx 0.87$, average fluid mass density $\rho_{b}=$351.3\pm
0.7\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$, and varying
Reynolds numbers.

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 $\hat{p}$ from the equation of state, Eq. 5. By
assuming $\hat{p}_{f}=p_{f}$, $\hat{p}=p_{f}$ and $\hat{\gamma}=\gamma$, the
individual contributions can be calculated. The outcome is illustrated in Fig.
9, again for a porosity $\phi\approx 0.87$, average fluid mass density
$\rho_{b}=$351.3\pm
0.7\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$. The Reynolds
number was now $\text{Re}=0.67\pm 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 $\gamma$
was calculated from the spherical mechanical pressure tensor. The value
increases monotonically with increasing density between $0.079\pm
0.009\text{\,}\mathrm{mN}\text{\,}{\mathrm{m}}^{-1}$ and $34\pm
1.\text{\,}\mathrm{mN}\text{\,}{\mathrm{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

$\hat{p}_{s}=p_{f}+\gamma A/V_{s}$ (24)

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 $\phi\nabla p_{f}=$-8629\pm
106.\text{\,}\mathrm{bar}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$$, in
integral solid pressure $(1-\phi)\nabla\hat{p}_{s}=$-9965\pm
215.\text{\,}\mathrm{bar}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$$, and
in surface tension $(A/V)\Delta\gamma=$370\pm
214.\text{\,}\mathrm{bar}\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{-1}$$. The
single parts sum to the total value shown in the figure. This value of
$\hat{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.

Figure 9: Profile of porosity $\phi\approx 0.87$, average fluid mass density
$\rho_{b}=$351.3\pm 0.7\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$
and Reynolds number $\text{Re}=0.67\pm 0.02$. Top: Visualization of the
simulation box. The border of a layer of volume $V_{l}$ is marked with a
dashed line. Bottom: Contributions from the bulk phases and the surface to the
integral pressure as function of the $x$-coordinate. The $x$-axis of the
visualization correspond to the graph. The points in the graph gives the
$x$-position of the center of the REV (the volume $V_{l}$).

A control of isothermal conditions was carried out. We know that a pressure
difference may generate a temperature difference, or vice versa [19]. 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 $\phi\approx[0.26,0.61]$ and the average bulk fluid mass
density was $\rho_{b}=$108\pm
1.\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$ in these plots.

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 10: The mass flux as a function of the negative integral pressure
gradient for porosities $\phi\approx[0.26,0.61]$ and average bulk fluid mass
density $\rho_{b}=$108\pm
1.\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}^{3}$$. The dashed lines
are linear fits to the points.

###  4.3 Conductivity and permeability

Fig. 8 gave 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 was plotted in Fig. 11 as a function of the factor
$\phi^{3}/(1-\phi)^{2}$ in the Kozeny-Carman equation. There was no convincing
relationship between the variables.

It is more interesting to consider the permeability $k$. Its value was
computed from the conductivity $l$ and the shear viscosity $\eta_{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 [35, 36]. 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}$.

Figure 11: Conductivity coefficient as a function of the porosity factor
$\phi^{3}/(1-\phi)^{3}$ for various average fluid mass densities. Figure 12:
Permeability as a function of inverse average pressure pressure.

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
$\tau=1$, shown as a grey dash-dotted line, does not bring out any new
physical insight.

Figure 13: Permeability in logarithmic scale as a function of the porosity
$\phi$ for varying average bulk fluid mass density. The Kozeny-Carman equation
is shown as a black dashed curve for a tortuosity $\tau=1$.

Fig. 14 shows the Klinkenberg coefficient $b$ as a function of porosity
$\phi$. 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.

Figure 14: Klinkenberg constant $b$ as a function of porosity $\phi$.

##  5 Conclusion and Perspectives

We have demonstrated a new thermodynamic procedure to find the local effective
pressure gradient that drives the mass flow through a porous medium. 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 as well as 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}$m2.

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 is the origin of Darcy scale behavior. For instance, we can define the
effective pressure without consideration of the capillary pressure. For two-
phase 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.

Competing Interests The authors have no relevant financial or non-financial
interests to disclose.

Author Contributions O.G. contributed to formal analysis, investigation,
methodology, software, visualization, conceptualization, writing original
drafts, reviewing, and editing. S.K. and D.B. contributed to formal analysis,
investigation, methodology, conceptualization, supervision, writing original
drafts, reviewing, and editing. T.T.T. computed viscosities and contributed to
reviewing and editing. M.S.B. and M.T.R. contributed to early versions of
software, visualization, methodology, conceptualization, reviewing, and
editing. All authors have read and agreed to the published version of the
manuscript.

Funding This work was funded by the Research Council of Norway through its
Centres of Excellence funding scheme, project number 262644, PoreLab.

Acknowledgments The authors acknowledge valuable comments to the last version
of the manuscript from Steffen Berg and Hamidreza Erfani Gahrooei. The
simulations were performed on resources provided by UNINETT Sigma2 - the
National Infrastructure for High-Performance Computing and Data Storage in
Norway with project numbers nn9229k and nn8022k. We thank the Research Council
of Norway for its Centres of Excellence funding scheme, project number 262644,
PoreLab.

Supplementary information to ”Local thermodynamic description of isothermal
single-phase flow in porous media”
Olav Galteland, Michael T. Rauter, Mina S. Bratvold, Thuat T. Trinh, Dick
Bedeaux and Signe Kjelstrup
PoreLab, Department of Chemistry, Norwegian University of Science and
Technology
March 4 2022

##  1 Surface tension

The fluid-solid surface tension was calculated from the spherical mechanical
pressure tensor of a single solid sphere surrounded by a fluid of varying mass
density. This was considered as the surface tension of a lattice with a large
lattice constant, such that the solid surfaces were far apart. The mechanical
pressure tensor can be written as a sum of an ideal gas contribution and a
virial contribution. The diagonal components of the mechanical pressure tensor
were calculated in spherical shells with origin at the center of the solid
sphere

$P_{\alpha\beta}(r)=\rho(r)k_{\text{B}}T\delta_{\alpha\beta}+P_{\alpha\beta}^{v}(r),$
(25)

where $r$ is the distance to the origin, $\rho(r)$ is the fluid mass density
of the spherical shell, $\delta_{\alpha\beta}$ is the Kroenecker delta and
$P^{s,v}_{\alpha\beta}$ is the virial contribution. The subscripts $\alpha$
and $\beta$ refer to the spherical coordinates, $(r,\theta,\phi)$. The virial
contribution is

$P_{\alpha\beta}^{v}(r)=-\frac{1}{V(r)}\sum_{i=1}^{N}\sum_{j>i}^{N}f_{ij,\alpha}\int_{C_{ij}\in
V(r)}\text{d}l_{\beta}.$ (26)

where $V(r)$ is the volume of the spherical shell, $f_{ij,\alpha}$ is the
$\alpha$-component of the force acting on particle $i$ due to particle $j$.
The line integral is along the Irving-Kirkwood contour $C_{ij}$ [Irving1950,
Harasima1958, Schofield1982], which is the straight line between particle $i$
and $j$. The line integral gives the length of the $\beta$-component of the
contour $C_{ij}$ that is contained in the spherical shell $V(r)$. The surface
tension is calculated from the pressure tensor

$\gamma=\frac{1}{R^{2}}\int_{R}^{r_{0}}(P_{N}-P_{T})r^{2}\text{d}r,$ (27)

where $P_{N}=P_{rr}$ is the normal to the fluid-solid surface and
$P_{T}=(P_{\theta\theta}+P_{\phi\phi})/2$ is tangential to it. The integral
limit $r_{0}$ is a position in the fluid far away from the fluid-solid surface
such that $P_{N}=P_{T}$. The surface tension is shown as a function of fluid
mass density in Fig. 15.

Figure 15: Solid-fluid surface tension as a function of fluid mass density for
a single solid sphere surrounded by a fluid.

## References

  * [1] T. Stylianopoulos, L. Munn, and R. Jain, “Reengineering the physical microenvironment of tumors to improve drug delivery and efficacy: From mathematical modeling to bench to bedside,” Trends in cancer, vol. 4, p. 292, 2018.
  * [2] A. Zlotorowicz, K. Jayasayee, P. Dahl, M. Thomassen, and S. Kjelstrup, “Tailored porosities of the cathode layer for improved polymer electrolyte fuel cell performance,” J. Power Sources, vol. 287, pp. 472–477, 2015\.
  * [3] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface.  Berlin: Springer, 1997.
  * [4] M. J. Blunt, Multiphase flow in permeable media: A pore-scale perspective.  Cambridge: Cambridge University Press, 2017.
  * [5] S. Whitaker, “Flow in porous media I: A theoretical derivation of Darcy’s law,” Transport in porous media, vol. 1, no. 1, pp. 3–25, 1986.
  * [6] K. Nordahl and P. S. Ringrose, “Identifying the representative elementary volume for permeability in heterolithic deposits using numerical rock models,” Mathematical geosciences, vol. 40, no. 7, pp. 753–771, 2008.
  * [7] W. Gray and M. Hassanizadeh, “Macroscale continuum mechanics for multiphase porous media flow including interface boundaries,” Adv. Water Resource, vol. 21, pp. 261–281, 1998.
  * [8] O. Galteland, D. Bedeaux, B. Hafskjold, and S. Kjelstrup, “Pressures inside a nano-porous medium. The case of a single phase fluid,” Front. Phys., vol. 7, p. 60, 2019.
  * [9] M. T. Rauter, O. Galteland, M. Erdős, O. A. Moultos, T. J. Vlugt, S. K. Schnell, D. Bedeaux, and S. Kjelstrup, “Two-Phase Equilibrium Conditions in Nanopores,” Nanomaterials, vol. 10, no. 4, p. 608, 2020.
  * [10] D. Bedeaux and S. Kjelstrup, “Fluctuation-dissipation theorems for flow in porous media,” Entropy, vol. 24, p. 46, 2021.
  * [11] S. Kjelstrup, D. Bedeaux, A. Hansen, B. Hafskjold, and O. Galteland, “Non-isothermal transport of multi-phase fluids in porous media. The entropy production,” Front. Phys., vol. 6, p. 126, 2018.
  * [12] S. Kjelstrup, D. Bedeaux, A. Hansen, B. Hafskjold, and O. Galteland, “Non-isothermal transport of multi-phase fluids in porous media. Constitutive equations,” Front. Phys., vol. 6, p. 150, 2019.
  * [13] B. D. Todd and P. J. Daivis, Nonequilibrium molecular dynamics: theory, algorithms and applications.  Cambridge: Cambridge University Press, 2017.
  * [14] B. Hafskjold, T. Ikeshoji, and S. K. Ratkje, “On the molecular mechanism of thermal diffusion in liquids,” Mol. Phys., vol. 80, no. 6, pp. 1389–1412, 1993.
  * [15] T. Ikeshoji and B. Hafskjold, “Non-equilibrium molecular dynamics calculation of heat conduction in liquid and through liquid-gas interface,” Mol. Phys., vol. 81, no. 2, pp. 251–261, 1994.
  * [16] B. Hafskjold and T. Ikeshoji, “Non equilibrium molecular dynamics simulation of coupled heat-and mass transport across a liquid/vapor interface,” Mol. Sim., vol. 16, no. 1-3, pp. 139–150, 1996.
  * [17] T. L. Hill, Thermodynamics of small systems, Parts I & II.  New York: Dover Publications, 1994.
  * [18] D. Bedeaux, S. Kjelstrup, and S. K. Schnell, Nanothermodynamics. General theory.  Trondheim: PoreLab, Norwegian University of Science and Technology, 2020\.
  * [19] M. T. Rauter, S. K. Schnell, B. Hafskjold, and S. Kjelstrup, “Thermo-osmotic pressure and resistance to mass transport in a vapor-gap membrane,” Phys. Chem. Chem. Phys., 2021.
  * [20] M. T. Rauter, S. K. Schnell, and S. Kjelstrup, “Cassie–Baxter and Wenzel States and the Effect of Interfaces on Transport Properties across Membranes,” J. Phys. Chem. B, 2021.
  * [21] M. Erdős, O. Galteland, D. Bedeaux, S. Kjelstrup, O. A. Moultos, and T. J. Vlugt, “Gibbs Ensemble Monte Carlo Simulation of Fluids in Confinement: Relation between the Differential and Integral Pressures,” Nanomaterials, vol. 10, no. 2, p. 293, 2020.
  * [22] O. Galteland, D. Bedeaux, and S. Kjelstrup, “Nanothermodynamic description and molecular simulation of a single-phase fluid in a slit pore,” Nanomaterials, vol. 11, no. 1, p. 165, 2021.
  * [23] O. Galteland, E. Bering, K. Kristiansen, D. Bedeaux, and S. Kjelstrup, “Legendre-Fenchel transforms capture layering transitions in porous media,” arXiv preprint arXiv:2111.15253, 2021.
  * [24] O. Galteland, M. T. Rauter, M. Bratvold, D. Bedeaux, and S. Kjelstrup, “Defining the pressures of a fluid in a nanoporous, heterogeneous medium,” arXiv preprint arXiv:2201.13060, 2022.
  * [25] J. Li, D. Liao, and S. Yip, “Coupling continuum to molecular-dynamics simulation: Reflecting particle method and the field estimator,” Physical Review E, vol. 57, no. 6, p. 7259, 1998.
  * [26] A. Stukowski, “Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool,” Model. Simul. Mater. Sc., vol. 18, no. 1, p. 015012, 2009.
  * [27] J. Kozeny, “Uber kapillare leitung der wasser in boden,” Sitzungsber. Akad. Wiss., vol. 136, pp. 271–306, 1927.
  * [28] P. C. Carman, “Fluid flow through granular beds,” T. I. Chem. Eng.-Lond., vol. 15, pp. 150–166, 1937.
  * [29] P. C. Carman, Flow of gases through porous media.  New York: Academic Press, 1956.
  * [30] C. F. Berg, “Permeability description by characteristic length, tortuosity, constriction and porosity,” Transport in porous media, vol. 103, no. 3, pp. 381–400, 2014.
  * [31] J. N. Israelachvili, Intermolecular and surface forces.  New York: Academic Press, 2015.
  * [32] R. C. Tolman, “The effect of droplet size on surface tension,” The journal of chemical physics, vol. 17, no. 3, pp. 333–337, 1949.
  * [33] A. Røsjorde, D. Fossmo, D. Bedeaux, S. Kjelstrup, and B. Hafskjold, “Nonequilibrium molecular dynamics simulations of steady-state heat and mass transport in condensation: I. local equilibrium,” Journal of colloid and interface science, vol. 232, no. 1, pp. 178–185, 2000.
  * [34] I. Inzoli, S. Kjelstrup, D. Bedeaux, and J.-M. Simon, “Transport coefficients of n-butane into and through the surface of silicalite-1 from non-equilibrium molecular dynamics study,” Microporous and Mesoporous materials, vol. 125, no. 1-2, pp. 112–125, 2009.
  * [35] L. Klinkenberg, “The permeability of porous media to liquids and gases,” American Petroleum Institute, pp. 200–213, 1941.
  * [36] A. Baehr and M. Hult, “Evaluation of unsaturated zone air permeability through pneumatic tests,” water Resources Research, vol. 27, no. 1, pp. 2605–2627, 1991.
  * [37] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, et al., “LAMMPS-A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Comput. Phys. Commun., p. 108171, 2021.
  * [38] W. A. Steele, The interaction of gases with solid surfaces.  Oxford: Pergamon Press, 1974.
  * [39] M. Miyahara and K. E. Gubbins, “Freezing/melting phenomena for Lennard-Jones methane in slit pores: a Monte Carlo study,” J. Chem. Phys., vol. 106, no. 7, pp. 2865–2880, 1997.
  * [40] A. McGaughey and M. Kaviany, “Thermal conductivity decomposition and analysis using molecular dynamics simulations. Part I. Lennard-Jones argon,” Int. J. Heat Mass Tran., vol. 47, no. 8-9, pp. 1783–1798, 2004.
  * [41] G. Rutkai, M. Thol, R. Span, and J. Vrabec, “How well does the Lennard-Jones potential represent the thermodynamic properties of noble gases?,” Mol. Phys., vol. 115, no. 9-12, pp. 1104–1121, 2017.
  * [42] B. Hafskjold, K. P. Travis, A. B. Hass, M. Hammer, A. Aasen, and Ø. Wilhelmsen, “Thermodynamic properties of the 3D Lennard-Jones/spline model,” Mol. Phys., vol. 117, no. 23-24, pp. 3754–3769, 2019.
  * [43] K. R. Kristiansen, “Transport Properties of the Simple Lennard-Jones/Spline Fluid I: Binary Scattering and High-Accuracy Low-Density Transport Coefficients,” Front. Phys., vol. 8, p. 271, 2020.
  * [44] S. H. Jamali, L. Wolf, T. M. Becker, M. de Groen, M. Ramdin, R. Hartkamp, A. Bardow, T. J. H. Vlugt, and O. A. Moultos, “Octp: A tool for on-the-fly calculation of transport properties of fluids with the order-n algorithm in lammps,” J. Chem. Inf. Model., vol. 59, pp. 1290–1294, 2019.
  * [45] W. Shinoda, M. Shiga, and M. Mikami, “Rapid estimation of elastic constants by molecular dynamics simulation under constant stress,” Phys. Rev. B, vol. 69, no. 13, p. 134103, 2004.
  * [46] D. Frenkel and B. Smit, Computational science series: Understanding molecular simulation: from algorithms to applications.  San Diego: Academic press, 2002.