Surface Permeability of Particulate Porous Media

Abstract

The dispersion process in particulate porous media at low saturation levels takes place over the surface elements of constituent particles and, as we have found previously by comparison with experiments, can be accurately described by superfast nonlinear diffusion partial differential equations. To enhance the predictive power of the mathematical model in practical applications, one requires the knowledge of the effective surface permeability of the particle-in-contact ensemble, which can be directly related with the macroscopic permeability of the particulate media. We have shown previously that permeability of a single particulate element can be accurately determined through the solution of the Laplace–Beltrami Dirichlet boundary value problem. Here, we demonstrate how that methodology can be applied to study permeability of a randomly packed ensemble of interconnected particles. Using surface finite element techniques, we examine numerical solutions to the Laplace–Beltrami problem set in the multiply-connected domains of interconnected particles. We are able to directly estimate tortuosity effects of the surface flows in the particle ensemble setting.

1 Introduction

Liquid transport in particulate porous media, such as sand, is customarily classified into fully saturated, funicular and pendular regimes of spreading (Bear 1972; Herminghaus 2005; Scheel et al. 2008a, b). The first two regimes of the liquid dispersion occur at relatively high saturation levels $s>s_c\approx 10\%$, where saturation s is defined as the ratio of the liquid volume $V_L$ to the volume of available voids $V_E$ in a sample volume element V, $s=\frac{V_L}{V_E}$. At high saturation levels, above the critical value $s_c$, liquid transport takes place in the pore space either fully or partially filled by the liquid.

Fig. 1

Illustration of the liquid distribution in particulate porous media (grey) with pendular rings (blue) at low saturation levels

Our prime concern here is the special case of liquid dispersion at low saturation levels. As the saturation level drops below the critical value, $s\le s_c$, that is to the value relevant to the pendular regime of spreading, the liquid volumes in the porous matrix become isolated (Herminghaus 2005; Scheel et al. 2008a, b). As a result, at low saturation levels, the liquid is only contained in the pendular rings formed at the locations of the particle contacts and on the particle rough surfaces, and the liquid transport can only occur over the matrix surface elements, as is illustrated in Fig. 1.

The analysis of this regime of wetting, which is crucial for studies of biological processes and spreading of nonvolatile liquids in arid natural environments and industrial installations, has shown that the liquid dispersion has many distinctive features and can be accurately described by the so-called superfast nonlinear diffusion equation (Lukyanov et al. 2012, 2019).

Our main concern here is the wetting cycle, when the liquid spreads over a dry porous matrix or over a matrix with a very low background saturation level up to $s_r\approx 2\%$. These conditions are similar to those in the case studied previously experimentally and theoretically in Lukyanov et al. (2019). The main driving force of the dispersion process, as is often the case during the wetting cycle, is capillary pressure developed at the moving front in the process of wetting of dry porous matrix, while the liquid bridges play a role of variable liquid reservoirs of uniform surface curvature. Note that, unlike at high saturation values $s>s_c$, at low saturation $s<s_c$ the capillary pressure is developed on the length scale of the surface roughness.

Theoretically, the superfast nonlinear diffusion equation belongs to a special class of mathematical models. Unlike in the standard porous medium equation (Vazquez 2006), where the coefficient of diffusion D(s) is either constant (leading to normal, Gaussian dispersive behaviour) or a monotonically vanishing function of saturation ($D(s)\propto s^m, \, m>0$, leading to hypodispersive behaviour), in this special case, the nonlinear coefficient of diffusion demonstrates divergent behaviour as a function of saturation, $D(s)\propto (s-s_0)^{-3/2}$, leading to hyperdispersive character of the spreading process. Here, $s_0$ is some minimal saturation level ($s_0\approx 0.5\%$), which could be only achieved in a state when the liquid bridges cease to exist completely (Lukyanov et al. 2012; Scheel et al. 2008a, b; Lukyanov et al. 2019). Note, in that respect, that in the domain of spreading liquid bridges are supposed to never vanish, so that the condition $s>s_0$ is always fulfilled in the model, and there is no actual singularity of the mathematical description (Lukyanov et al. 2012, 2019). The character of the dispersion process may often change from hypodispersive to hyperdispersive depending on the spreading scenario, that is on the driving forces of the spreading process (de Gennes 1985; Bacri et al. 1985; Novy et al. 1989; Toledo et al. 1993; Meleá 2003), especially in multi-phase systems during forced imbibition cycles used in petroleum engineering (Meleá 2003). If the forcing, external pressure at the system inlet is strong enough, the low saturation zone may not have a chance to be developed, and hypodispersive character of the spreading front is observed.

Returning to our mathematical model at low saturation values, in the macroscopic approximation, that is after averaging over some volume element containing many particles of the porous medium, the diffusion process in the slow creeping flow conditions can be described by the following nonlinear diffusion equation

$$\begin{aligned} \frac{\partial s}{\partial t}= \nabla \cdot \left\{ D(s) \nabla s \right\} , \quad t>0, \end{aligned}$$

(1)

where

$$\begin{aligned} D(s)= \frac{D_0(s)}{(s-s_0)^{3/2}}, \quad s>s_0, \end{aligned}$$

for $D_0 > 0$.

The details of derivation of (1) can be found in Lukyanov et al. (2012) and Lukyanov et al. (2019), here we note that, the resultant governing nonlinear Eq. (1) directly follows from the conservation of mass principle

$$\begin{aligned} \frac{\partial (\phi s)}{\partial t} + \nabla \cdot \mathbf{Q}=0, \end{aligned}$$

(2)

where $\phi$ is porosity defined as $\displaystyle \phi =\frac{V_E}{V}$, which is further assumed to be constant, and $\mathbf Q$ is the macroscopic flux density. The macroscopic flux density $\mathbf Q$ is defined in such a way that the total flux through the surface of a macroscopic sample volume element is given by the surface integral $\int \mathbf{\!\! Q\cdot n}\, \hbox {d}S$, where $\mathbf n$ is the normal vector to the surface of the sample volume element.

To obtain (1) from (2), one needs to apply the capillary pressure–saturation relationship (Halsey and Levine 1998; Lukyanov et al. 2012, 2019) dictated by the liquid bridges behaviour

$$\begin{aligned} p = -p_0\frac{A_c}{(s-s_0)^{1/2}} \end{aligned}$$

(3)

and the local Darcy’s law (Rye et al. 1998; Or and Tuller 2000) describing the surface flow in the rough layer of the particle elements

$$\begin{aligned} -\frac{\kappa _m}{\mu } \nabla u =\mathbf{q}. \end{aligned}$$

(4)

here $A_c=\sqrt{\frac{3}{4}\, \frac{1-\phi }{\phi } \frac{N_c}{\pi }}$, $N_c$ is the coordination number, that is the average number of bridges per a particle, $p_0=\frac{2\gamma }{R}\cos \phi _c$, $\gamma$ is the coefficient of the surface tension of the liquid, $\phi _c$ is the contact angle made by the free surface of the liquid bridge with the rough solid surface of the constituent particles, R is an average radius of the porous medium particles, $\mathbf{q}$ and u are the averaged local flux density and pressure in the rough surface layer, $\mu$ is liquid viscosity, and $\kappa _m$ is the local coefficient of permeability of the rough surface, which is proportional to the average amplitude of the surface roughness $\delta _R$, that is the width of the surface layer conducting the liquid flux

$$\begin{aligned} \kappa _m\propto \delta _R^2. \end{aligned}$$

(5)

One needs to emphasise here that two levels of averaging are involved in obtaining the final governing Eq. (1). While Eqs. (1)–(3) are ‘truly’ macroscopic, that is obtained by averaging using a volume element V containing many grain particles, Eq. (4) is only an average over some rough area of a single particle containing many surface irregularities, so that quantities $\mathbf{q}$ and u are also only local averages over that sample surface area.

Therefore, to transit from (4) to the macroscopic description, the spatial averaging theorem formulated in Whitaker (1969) should be applied. That is, using intrinsic liquid averaging $\langle \cdots \rangle ^l=V_l^{-1}\!\int _{V_l} \ldots \hbox {d}^3 x$, where $V_l$ is liquid volume within the sample volume V, one has $\langle u \rangle ^l=p$ and $\langle \mathbf{q}\rangle ^l \frac{S_e}{S} = \mathbf{Q}$. Here, S is the surface area of the sample volume V with the effective area of entrances and exits $S_e$. Note, the ratio $S_e/S$ is not just a geometric property, but also takes into account the connectivity of the porous elements. For example, the effective area of entrances and exits $S_e$ is only defined by the pathways open to the flow.

As a result of the two-level averaging

$$\begin{aligned} D_0(s)=A_c \frac{K(s)}{\mu }\frac{p_0}{2\phi }, \end{aligned}$$

where $K(s)=\kappa _m \frac{S_e}{S}$ is the coefficient of permeability defined by

$$\begin{aligned} \mathbf{Q} = -\frac{K}{\mu }\nabla p. \end{aligned}$$

Fig. 2

Illustration of the flow and solution domains on the surface $\varGamma$ of a spherical particle, and their geometric arrangements. In the picture, $\varOmega _0$ is the domain of the surface flow and the surface area covered by the liquid bridges corresponds to the domains $\varOmega _1$ and $\varOmega _2$

The global surface permeability of the particles K as a function of saturation is one of the main elements of the model to accurately represent liquid dispersion at low saturation levels. It is fully defined by the particle geometry and the geometry of the liquid bridge contact areas, Figs. 1 and 2.

In particular, the disposition and the size of the liquid bridges on the particle surface, that is the size of the domains $\varOmega _{1,2}$ and the angle $\alpha$, should play a leading role in defining the resistance to the surface flow. It is not difficult to discern that any variations of the contact area covered by the liquid bridges (pendular rings), that is areas $\varOmega _{1,2}$ shown in Fig. 2, or the value of the bridge volume, should affect the global permeability.

Previously, we have shown that permeability of a single particle element can be determined by means of a solution to the equivalent Laplace–Beltrami boundary value problem formulated in the flow domain $\varOmega _0$ with the boundaries $\partial \varGamma _{1,2}$ in Fig. 2 (Sirimark et al. 2018). We briefly formulate that problem and summarise the previous results in the next part. Here, we note that, based on the analysis of the problem, we have been able to show that in a special azimuthally symmetric case of spherical particles, when the two areas covered by the liquid bridges, domains $\varOmega _1$ and $\varOmega _2$ in Fig. 2, are oriented symmetrically to each other, that is at $\alpha =\pi$, the permeability K is supposed to follow the scaling

$$\begin{aligned} K(s)\propto \frac{1}{|\ln (s-s_0)|}. \end{aligned}$$

We have studied several generalisations of the symmetric problem, such as arbitrary oriented domains, $\alpha \ne \pi$, on the surface of the spherical particle, and particles of arbitrary shapes emulating the shape of a real sand grain. While variations of the particle shape were found to produce a relatively modest effect on the particle surface permeability, the orientation of the boundaries, emulating tortuosity effects, was found to produce a stronger impact due to the substantial variation of the distance, on average, between the boundary contours $\partial \varGamma _{1,2}$. It became clear that while the previously obtained scaling was a good first step to estimate the surface permeability of particulate porous media, a more general case of an ensemble of interconnected particles should be analysed to enhance the model predictive power and at the same time to estimate rigorously the effects of tortuosity of the surface flow in the particle assembly. In this study, we will simulate a general case of an ensemble of many particles linked by liquid bridges. We will concentrate on the bunch of spherical particles, but of different radii and randomly arranged in configurations. We compare the random pack configuration results with some symmetric case to estimate the effects of tortuosity and formulate practical recipes to apply the superfast diffusion model.

2 Microscopic Model of the Surface Permeability of the Elements

Microscopically, the liquid creeping flow through the surface roughness of each particle can be described by a local Darcy-like relationship (4) between the surface flux density $\mathbf{q}$ and averaged (over some area containing many surface irregularities) pressure in the grooves u (Rye et al. 1998; Or and Tuller 2000). Assuming incompressibility of the liquid and that the liquid layer thickness is constant $\delta _R=\mathrm{const}$, one has

$$\begin{aligned} \nabla \cdot \mathbf{q}=0. \end{aligned}$$

(6)

Equation (4) taking into account (6) can then be transformed into the Laplace–Beltrami equation defined on the surface $\varGamma$ of the particle

$$\begin{aligned} \Delta _{\varGamma } u =0. \end{aligned}$$

(7)

here $\Delta _{\varGamma }$ designates the Laplace–Beltrami operator, which is defined on the surface element $\varGamma$ through the surface gradient $\nabla _{\varGamma }$ tangential to the surface. Formally, let $\mathbf{{n}}_{\varGamma }$ denote the unit normal to the surface $\varGamma$, Fig. 2. Then, one can define the surface gradient of a smooth function u as $\nabla _{\varGamma } u:=\nabla u - (\nabla u \cdot \mathbf{{n}}_{\varGamma }) \mathbf{{n}}_{\varGamma }$, and then, the Laplace–Beltrami operator is defined as $\Delta _{\varGamma } u = \nabla _{\varGamma }\cdot \nabla _{\varGamma } u$.

The second assumption $\delta _R=\hbox {const}$ implies that the surface layer is fully saturated, that is its content is not changing on the particle surface. The approximation of the fully saturated rough surface layer is well fulfilled, if the characteristic pressure amplitude |u| is less than the capillary pressure amplitude defined on the length scale of the surface roughness $\delta _R$, which is of the order of $\delta _R\sim 1\,\upmu \text{ m }$ in typical sands (Alshibli and Alsaleh 2004), as is demonstrated in Rye et al. (1998). That is, $|u| < u_c=\frac{\gamma }{\delta _R}$, and, for example for water ($\gamma =72\, \text{ mN }/\text{ m }$) at $\delta _R =1\,\mu \text{ m }$, this results in $|u| < 7.2\times 10^4\, \text{ Pa }$.

Alternatively, if the surface layer somehow is not fully saturated, parameter $\delta _R$ should be interpreted as the characteristic width of the liquid layer within the rough surface layer and one needs to presume that variations of the pressure $|\delta u|$ are negligible $|\delta u|\ll u_c$. This is usually the case in slow, creeping flow conditions in porous media, and in fact, it is a criterion for the use of macroscopic approximation to such flows (Bear 1972). As is shown in Lukyanov et al. (2019), strong negative capillary pressure on the level of $u_c$ is only expected at the moving front, so that the approximation is well fulfilled in the macroscopic flow domain. Note also that, it is always assumed throughout this study that

$$\begin{aligned} \delta _R\ll R, \end{aligned}$$

that is the amplitude of the surface roughness (or the width of the liquid layer) is always much smaller than the particle size.

2.1 Permeability of a Single Particle Element

Consider, as the simplest example, a spherical particle of radius R with a closed surface $\varGamma$, which is split into three sub-domains $\varOmega _0$, $\varOmega _1$ and $\varOmega _2$ with the surface boundaries between them $\partial \varGamma _1$ and $\partial \varGamma _2$, as is shown in Fig. 2. The location of the sub-domains $\varOmega _1$ and $\varOmega _2$ to each other on the surface is fixed by the tilt angle $\alpha$. The sub-domains $\varOmega _1$ and $\varOmega _2$ correspond to the contact area covered by the liquid in the bridges, while the surface flow, described by (4), takes place in $\varOmega _0$. Our prime concern is permeability of the surface elements, so that we only consider steady-state problems.

The distribution of liquid pressure u, as it follows from (7), should satisfy the Laplace–Beltrami equation now defined on the surface of the sub-domain $\varOmega _0$

$$\begin{aligned} \Delta _{\varOmega _0} u =0. \end{aligned}$$

(8)

Note that, in fact, the condition of the fully saturated surface layer is not essential in calculation of the flows over one particle element of the porous media. It is sufficient to presume that the variation of the capillary pressure on the length scale of the particle $|\delta u|$ is negligible, that is $|\delta u | \ll u_c$. In the case, when the surface layer is not fully saturated, parameter $\delta _R$ should be interpreted as the effective thickness of the layer filled by the liquid.

At the same time, liquid pressure variation in the bridges is negligible in slow, creeping flows in comparison with that in $\varOmega _0$. So that, one can assume that

$$\begin{aligned} \left. u \right| _{\partial \varGamma _1}=U_1=\hbox {const}, \quad \left. u\right| _{\partial \varGamma _2}=U_2=\hbox {const}, \end{aligned}$$

(9)

which are the boundary conditions to the Laplace–Beltrami Dirichlet boundary value problem. The Dirichlet boundary value problem (8)–(9) has at least a unique weak solution, if the domain $\varOmega _0$ and the boundaries $\partial \varGamma _{1,2}$ are smooth (Dziuk 1988; Pigola et al. 2005; Dziuk and Elliott 2013), which, if it is found, allows to calculate the total flux through the particle element

$$\begin{aligned} Q_\mathrm{T}=\delta _R\frac{\kappa _m}{\mu }\int _{\partial \varGamma _1} \frac{\partial u}{\partial n_s}\, \mathrm{d}l = -\delta _R\frac{\kappa _m}{\mu } \int _{\partial \varGamma _2} \frac{\partial u}{\partial n_s} \, \mathrm{d}l, \end{aligned}$$

(10)

where $\mathbf n_s$ is the normal vector to the domain boundaries $\partial \varGamma _{1,2}$ on the surface, $\delta _R$ is the average amplitude of the surface roughness, that is the width of the surface layer conducting the liquid flux and the line integral is taken along a closed curve in $\varOmega _0$, for example the boundary $\partial \varGamma _1$.

If the total flux $Q_\mathrm{T}$ is determined, one can define the global permeability coefficient of a single particle $K_1$. This can be done, if we assume that the particle has a characteristic size D and so that it can be enclosed in a volume element $V=D^3$ with the characteristic side surface area $D^2$. Then, the effective flux density Q can be represented in terms of $K_1$ (and the total flux $Q_\mathrm{T}$)

$$\begin{aligned} Q=\frac{Q_\mathrm{T}}{D^2}=-\frac{K_1}{\mu }\frac{U_2-U_1}{D}, \end{aligned}$$

(11)

if the flow is driven by the constant pressure difference $U_2-U_1$ applied to the sides of the volume element.

2.2 Surface Permeability of a Sphere in the Case of Azimuthally Symmetric Domain Boundaries

Consider now a spherical particle in an azimuthally symmetric case, when the domain boundaries $\partial \varGamma _1$ and $\partial \varGamma _2$ are oriented at the reflex angle $\alpha =\pi$ and have a circular shape. We use a spherical coordinate system with its origin at the particle centre and the polar angle $\theta$ counted from the axis of symmetry passing through the centre of the circular contour $\partial \varGamma _1$. In this case, the Dirichlet boundary value problem (8)–(9) admits an analytical solution, so that particle permeability can be determined explicitly. Indeed, problem (8)–(9), if we assume that the liquid pressure distribution u is a function of $\theta$ only and independent of the azimuthal angle, is equivalent to

$$\begin{aligned} \frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left( \sin \theta \frac{\partial u}{\partial \theta }\right) =0, \quad \theta _0< \theta < \pi -\theta _1, \end{aligned}$$

(12)

with the boundary conditions

$$\begin{aligned} \left. u\right| _{\theta =\theta _0}=U_1, \quad \left. u \right| _{\theta =\pi -\theta _1}=U_2. \end{aligned}$$

(13)

The analytic solution to problem (12)–(13) after applying the boundary conditions can be represented in the following form

$$\begin{aligned} u=\varPsi _0 (U_2-U_1) \ln \left\{ \frac{\sin \theta }{\sin \theta _0} \frac{1+\cos \theta _0}{1+\cos \theta } \right\} +U_1, \end{aligned}$$

(14)

where

$$\begin{aligned} \varPsi _0=\frac{1}{\ln \left\{ \frac{\sin \theta _1}{\sin \theta _0} \frac{1+\cos \theta _0}{1-\cos \theta _1} \right\} }. \end{aligned}$$

One can now calculate the total flux and the permeability, using its definition (11),

$$\begin{aligned} Q_\mathrm{T}=-\frac{K_1}{\mu } D (U_2-U_1)=-2\pi \sin \theta _0 \delta _R \frac{\kappa _m}{\mu }\left. \frac{\partial u}{\partial \theta }\right| _{\theta =\theta _0} \end{aligned}$$

$$\begin{aligned} =-(U_2-U_1) 2\pi \delta _R \varPsi _0 \frac{\kappa _m}{\mu }. \end{aligned}$$

(15)

So that, taking $D=2R$,

$$\begin{aligned} K_1=\pi \varPsi _0 \frac{\delta _R}{R} \kappa _m. \end{aligned}$$

(16)

Parametrically, the coefficient of permeability (16) is inversely proportional to the particle radius R, so that larger particles create stronger resistance to the flow. Noticeably, the coefficient demonstrates strong dependence on the surface layer thickness $\delta _R$, that is $K_1\propto \delta _R^3$ since it is anticipated that $\kappa _m\propto \delta _R^2$, so that evaluation of this parameter in applications is crucial for the accurate estimates of the liquid dispersion rates.

One can see, if we take $\theta _1=\theta _0$, in fact assuming small variations of the bridge size and the pressure over one particle diameter, and $\theta _0\ll 1$, in fact considering small values of saturation, $s\ll 1$, that the permeability coefficient $K_1$ tends to zero, that is

$$\begin{aligned} K_1 = \frac{\delta _R}{2R} \frac{\pi \kappa _m}{|\ln \theta _0 |} + o\left( \frac{1}{|\ln \theta _0|}\right) . \end{aligned}$$

(17)

How does the result affect the superfast diffusion model (1), and basically how can it be incorporated into the main diffusion equation? If we approximate the permeability coefficient K by $K_1$ obtained in the azimuthally symmetric case at $\theta _1=\theta _0$, (17), and, using an approximate relationship between the radius of curvature $R\sin \theta _0$ of the boundary contour $\partial \varGamma _1$ and the pendular ring volume (Herminghaus 2005), one can show that

$$\begin{aligned} \sin ^2\theta _0 \approx \theta _0^2 = \sqrt{s-s_0}. \end{aligned}$$

Therefore, finally

$$\begin{aligned} K(s)\approx K_1(s)\approx 2 \frac{\delta _R}{R} \frac{\pi \kappa _m}{|\ln (s-s_0)|}. \end{aligned}$$

(18)

As it follows from (18), the distinctive particle shape results in logarithmic correction to the main nonlinear superfast diffusion coefficient $D(s)=\frac{D_0(s)}{(s-s_0)^{3/2}}$, such that

$$\begin{aligned}D(s)\propto \frac{1}{|\ln (s-s_0)|(s-s_0)^{3/2}}. \end{aligned}$$

Apparently, the correction will mitigate to some extent the divergent nature of the dispersion at the very small saturation levels $s\approx s_0$, smoothing out the characteristic dispersion curves.

Fig. 3

Illustration of the solution domains in a system of two coupled spherical particles and their geometric arrangements

2.3 Surface Permeability of a Chain of Spheres in the Case of Azimuthally Symmetric Domain Boundaries

Consider now how the problem can be formulated in the case of several particles arranged in a single chain, as is illustrated in Fig. 3 in the case of two coupled by the bridge particles. To create the flow in the system of two coupled particles, one can set pressure difference between $\partial \varGamma _1^{(1)}$ and $\partial \varGamma _2^{(2)}$. Mathematically, this is equivalent of setting Dirichlet boundary conditions on $\partial \varGamma _1^{(1)}$ and $\partial \varGamma _2^{(2)}$ as in the previous case of a single particle. The boundaries $\partial \varGamma _1^{(2)}$ and $\partial \varGamma _2^{(1)}$ are ‘internal’, that is common to the bridge linking the flow between the two particles. Apparently, the pressure is supposed to be the same on the two contours

$$\begin{aligned} u_1\left. \right| _{\partial \varGamma _1^{(2)}}=u_2\left. \right| _{\partial \varGamma _2^{(1)}}=\mathrm{const} \end{aligned}$$

(19)

and due to conservation of mass in steady-state conditions in the absence of sinks and sources of the liquid one has

$$\begin{aligned} \oint _{\partial \varGamma _1^{(2)}} \nabla u_1 \cdot \mathbf{n}_{s_1}|_{\partial \varGamma _{1}^{(2)}}\, \mathrm{d}l =-\oint _{\partial \varGamma _2^{(1)} }\nabla u_2 \cdot \mathbf{n}_{s_2}|_{\partial \varGamma _{2}^{(1)}}\, \mathrm{d}l \end{aligned}$$

(20)

where $\mathbf{n}_{s_1}$ and $\mathbf{n}_{s_2}$ are the outward tangential normal vectors to the boundary contours $\partial \varGamma _{1,2}^{(2,1)}$, and $u_1$ and $u_2$ designate distribution of pressure on each particle, respectively.

As a result, the problem to define the flow and the permeability of the system corresponds to a system of two Laplace–Beltrami equations

$$\begin{aligned} \frac{1}{\sin \theta } \frac{\partial }{\partial \theta } \left( \sin \theta \frac{\partial u_1}{\partial \theta } \right) = 0, \quad \theta _0 \le \theta \le \pi - \theta _0 \end{aligned}$$

(21)

and

$$\begin{aligned} \frac{1}{\sin \theta } \frac{\partial }{\partial \theta } \left( \sin \theta \frac{\partial u_2}{\partial \theta } \right) = 0, \quad \theta _0 \le \theta \le \pi - \theta _0, \end{aligned}$$

(22)

but with a slightly different set of the boundary conditions

$$\begin{aligned} u_1 \left. \right| _{\theta =\theta _0} = U_1 \end{aligned}$$

(23)

$$\begin{aligned} u_2 \left. \right| _{\theta =\pi -\theta _0} = U_2, \end{aligned}$$

(24)

$$\begin{aligned} u_1 \left. \right| _{\theta =\pi -\theta _0} = u_2 \big |_{\theta =\theta _0} \end{aligned}$$

(25)

and

$$\begin{aligned} \left. \left( \sin \theta \dfrac{\partial u_1}{\partial \theta } \right) \right| _{\theta =\pi -\theta _0} = \left. \left( \sin \theta \dfrac{\partial u_2}{\partial \theta }\right) \right| _{\theta =\theta _0}, \end{aligned}$$

(26)

where $\theta _0$, as before, defines the size of the bridge footprint on the particle surface in the spherical coordinate system with the axis of symmetry passing through the centre of the bridge area, Fig. 3. Since we assumed, due to relatively small variations of pressure over a few grain particles, that all bridges are roughly identical, we have only one parameter $\theta _0$ to describe the bridge size.

Apparently, Eqs. (21) and (22) can be integrated twice, similar to the previous problem of a single particle (12), to obtain

$$\begin{aligned} u_{1} = C_0 \ln \frac{\sin \theta }{1+\cos \theta } + C_1, \end{aligned}$$

(27)

$$\begin{aligned} u_{2} = B_0 \ln \frac{\sin \theta }{1+\cos \theta } + B_1, \end{aligned}$$

(28)

where $C_{1,2}$ and $B_{1,2}$ are free constant parameters to be found from the boundary conditions.

It is not difficult to see from (26), that one has $C_0=B_0$ implying continuity of the contact flux. Applying the remaining boundary conditions (23)–(25), from (27) and (28)

$$\begin{aligned} u_{1} = \varPsi _0^{(2)} (U_2-U_1) \ln \left( \frac{\sin \theta }{1+\cos \theta } \frac{1+\cos \theta _0}{\sin \theta _0} \right) + U_1, \end{aligned}$$

(29)

$$\begin{aligned} u_{2} = \varPsi _0^{(2)} (U_2-U_1) \ln \left( \frac{\sin \theta }{1+\cos \theta } \frac{1 - \cos \theta _0}{\sin \theta _0} \right) + U_2, \end{aligned}$$

(30)

where

$$\begin{aligned} \varPsi _0^{(2)}=\frac{1}{2\ln \left( \frac{1+\cos \theta _0}{1-\cos \theta _0}\right) }. \end{aligned}$$

One can now calculate total flux and define permeability of the coupled spherical particles $K_2$

$$\begin{aligned} Q_\mathrm{T}= & {} -\frac{K_2}{2\mu } D (U_2-U_1)=-2\pi \sin \theta _0 \delta _R \frac{\kappa _m}{\mu }\left. \frac{\partial u_1}{\partial \theta }\right| _{\theta =\theta _0}\nonumber \\= & {} -(U_2-U_1) 2\pi \delta _R \varPsi _0^{(2)} \frac{\kappa _m}{\mu }, \end{aligned}$$

(31)

where D is the characteristic length scale of the cross section in the problem, $\kappa _m$ is local permeability of the surface layer, and $\delta _R$ is the layer width and $\mu$ is liquid viscosity.

So that, taking simply $D=2R$,

$$\begin{aligned} K_2=2\pi \varPsi _0^{(2)} \frac{\delta _R}{R} \kappa _m. \end{aligned}$$

(32)

One can see that the permeability of a system of two coupled particles $K_2$ is identical to that of a single particle (16), basically from Eqs. (16) and (32)

$$\begin{aligned} \frac{K_2}{K_1}=\frac{2\varPsi _0^{(2)}}{\varPsi _0}=1. \end{aligned}$$

It is not difficult to discern by deduction that in a general case of N coupled particles in a chain

$$\begin{aligned} K_N=\pi N \varPsi _0^{(N)} \frac{\delta _R}{R} \kappa _m=K_1 \end{aligned}$$

(33)

where

$$\begin{aligned} \varPsi _0^{(N)}=\frac{1}{N\ln \left( \frac{1+\cos \theta _0}{1-\cos \theta _0}\right) }. \end{aligned}$$

Note, a setup of many beads coupled by liquid bridges can be potentially used in microfluidics to create flexible water channels (Chen et al. 2016). If the radius of curvature of the particle chain is much larger than the particle size, the transport through such a microfluidic system should be defined by the permeability of a single particle, relationship (16), if the particle shape can be approximated by a sphere.

Fig. 4

Schematic illustration of the particle ensemble and the sample volume element setup

One can conclude in this part, that if the porous media configuration is made of parallel chains of particles oriented symmetrically to each other, and the flow is generated along the chains, the surface permeability given by Eq. (16) is the exact result.

3 Surface Permeability of a Randomly Packed Particle Ensemble

In real systems, the particles are interconnected randomly, so that the effects of tortuosity should substantially affect the permeability of the system (Carman 1937; Lorenz 1961; Bear 1972; Ghanbarian et al. 2013). To analyse those effects, we consider an ensemble of spherical particles randomly packed, as is shown in Fig. 4. The randomly packed configuration of approximately 3000–7000 particles has been generated by means of a molecular dynamics technique by applying a constant force to every particle placed in a box with reflecting boundaries (in the perpendicular direction to the box side), and interacting via the Lennard–Jones potential with different characteristic length scales R distributed normally, that is with the probability of the particle radius $W(R)\propto \exp \left( -\frac{(R-R_0)^2}{\Delta R^2}\right)$ at $\Delta R/R_0 = 0.3$. In this study, there were particles with three different characteristic dimensions $R_1=1.3\, R_0$, $R_2=R_0$ and $R_3=0.7\, R_0$. The resultant porosity in the configurations was about $48\%$.

Table 1 Tabulated values of the system parameters in the random configurations used in the study

To obtain the configuration, the particle temperature controlled by the thermostat has been gradually reduced to bring the system to a minimum energy, frozen state. A representative sample volume element with dimensions $L_x^B, L_y^B, L_z^B$ then was cut off the system, as is illustrated in Fig. 4, containing $N_S=13{-}17$ particles, see Table 1 for details. We have generated several statistically independent sample configurations, and, as in the previous examples, set constant pressure difference $U_2-U_1$ at the boundaries of the sample elements, Figs. 4 and 5.

The Laplace–Beltrami method then has been applied after establishing the position of the liquid bridges coupling the particles in the sample. Two particles (of radii $R_1$ and $R_2$) are assumed to be coupled by a liquid bridge if the distance between their centres r was only slightly larger than the sum of their radii

$$\begin{aligned} R_1+R_2 \le r < R_1 + R_2 +0.05\max (R_1,R_2). \end{aligned}$$

Fig. 5

Illustration of the particle sample and the flow domains

Fig. 6

Illustration of the pendular ring characteristic geometry

The size of a single liquid bridge footprint $H_B$ on the particle surface can be characterised, as before, by the polar angle $\theta _0$ in the polar coordinate system with the symmetry axis passing through the centre of the circular contour, the boundary of the area covered by the bridge, as is shown in the symmetric case in Fig. 3. That is, $H_B=2R_k\sin \theta _0^{(k)}$. Due to the specific geometric properties of the pendular rings (constant mean curvature surface), we assume that even in the case of a distribution of particles with different radii $R_k$, the size of the bridge area in the sample is approximately the same in the low saturation limit $s\ll 1$ ($\theta _0^{(k)}\ll 1$) (Herminghaus 2005; Halsey and Levine 1998; Willett et al. 2000).

Indeed, when $s\ll 1$, the pressure in the pendular ring p is defined by the smallest radius of curvature $r_1$, Fig. 6, $p\approx - \gamma \cos \phi _c/r_1$, which is related with the second radius

$$\begin{aligned} r_1\approx r_2^2/2R_k, \end{aligned}$$

(34)

so that when $s\ll 1$, one has $r_2\gg r_1$. Obviously, $r_2$ defines the size of the area covered by the bridge, $H_B=2r_2=2R_k\sin \theta _0^{(k)}$.

If we have two particles of different radii, say $R_1$ and $R_2$, in contact, the size of the bridge area will be approximately the same $r_2^{(1)}\approx r_2^{(2)}$ at low saturation levels, $s\ll 1$, with the difference being proportional to $r_1$, that is

$$\begin{aligned} \frac{ r_2^{(1)} - r_2^{(2)} }{ \max (R_1,R_2) }=O\left( \frac{r_1}{\max (R_1,R_2)} \right) . \end{aligned}$$

(35)

Apparently, in a general case, no analytic solution is expected to the Laplace–Beltrami problem and a well established surface finite element technique (Dziuk 1988; Pigola et al. 2005; Dziuk and Elliott 2013) is applied after the tessellation of the domains, as is shown in Fig. 7 for one particulate element with two boundary contours. The numerical method has been validated against analytical solutions previously demonstrating prescribed order of accuracy and numerical convergence, see details in Sirimark et al. (2018).

The number of particles in the sample volume element was negotiated between computational efficiency of the surface finite element method (so that, practically, any mesh resolution can be used to deal with any details on the boundary contours $\partial \varGamma _{k}^{(l)}$ and on the particle surfaces) and fluctuations of the averaged quantities obtained using the sample element, which are proportional to $N_S^{-1/2}\approx 25\%$. Moderate increase of the number of particles in the sample may significantly increase computational time to obtain highly resolved numerical solutions, while at the same time would not substantially reduce the effect of particle number fluctuations.

As one can see, problems (21)–(26) and hence total flux $Q_\mathrm{T}$ through a particle or a chain of particles, (15) or (31), are invariant under the transformation of the particle dimension R provided that the angular size of the bridge $\theta _0$ is fixed. In what follows, we change to non-dimensional description by normalising length scales by the average radius $R_0$ of the particles in the sample and pressure by the characteristic capillary pressure $p_0=2\gamma \cos \phi _c/R_0$. The flux $Q_\mathrm{T}$ will be normalised by the characteristic value

$$\begin{aligned} Q_0=p_0 \delta _R \frac{\kappa _m}{\mu }\frac{\bar{U}_2-\bar{U}_1}{\bar{L}^B_z} \end{aligned}$$

inspired by the analytical result (15) and by the non-dimensional sample box surface area $S_0=\bar{L}_x^B \bar{L}_y^B$, where non-dimensional quantities $\bar{L}_{x,y,z}^B=L_{x,y,z}^B/R_0$ and $\bar{U}_{1,2}=U_{1,2}/p_0$. The latter normalisation allows to bring simulation results in slightly different geometric settings, as is detailed in Table 1, into equivalent conditions suitable for comparison, that is basically providing the non-dimensional permeability $\bar{K}=\frac{K}{\kappa _m}\frac{R_0}{\delta _R}$.

Fig. 7

Illustration of the tessellated flow domain of a particle for the surface finite element method

Fig. 8

Reduced total flux $Q_\mathrm{T}/S_0 Q_0$, $\displaystyle Q_0=p_0 \delta _R \frac{\kappa _m}{\mu }\frac{\bar{U}_2-\bar{U}_1}{\bar{L}^B_z}$ and $S_0=\bar{L}_x^B \bar{L}_y^B$, as a function of $\displaystyle \varPsi _0(H_B/2R_0)$ as is defined in (38). The error bar indicates the statistical error, which is expected due to the fluctuations of the number of particles in the samples

Schematically, the simulation domains for the Laplace–Beltrami problem are shown in Fig. 5. As in the previous case of particles coupled in a chain, there are internal, common boundaries, contours $\partial \varGamma _{k}^{(l)}$, $k\ne l$, where the continuity boundary conditions are applied and external boundaries, contours $\partial \varGamma _k^{(k)}$, where the Dirichlet boundary conditions are set to generate a flow through the system. The pressure value on the contours facing the bottom of the simulation box (for example, $\partial \varGamma _2^{(2)}, \partial \varGamma _3^{(3)}$ and $\partial \varGamma _4^{(4)}$ in Fig. 5) is set to $U_2$ and on the contours facing the top side of the pack (for example, $\partial \varGamma _0^{(0)}$ and $\partial \varGamma _1^{(1)}$ in Fig. 5) is set to $U_1$. The values of the boundary pressure $U_{1,2}$ were identical in the simulations involving different configurations.

Geometrically, the external boundary contours are oriented in the flow direction, as is illustrated in Fig. 5. While this particular orientation seems to be arbitrary or may even look artificial, within the statistical approach, the choice of the boundary contour orientation should not render any excessive (in excess of the statistical errors due to the particle number fluctuations) influence upon the results, that is the value of the total flux and the ‘macroscopic’ permeability. A posteriori, one can see that this seemed to be the case, Fig. 8, as in different configurations, Table 1, the resultant curves are close and parallel to each other.

There are two main questions, we would like to answer in this part of the study. First, how does permeability of the particle sample depend on the composition? Basically, how strong are there fluctuations? Secondly, what is the contribution of the tortuosity effects? To obtain statistically meaningful results, we consider several randomly generated configurations, as is summarised in Table 1. We would like to stress here, that all configurations have been cut off from statistically independent particle distributions generated with the help of random initial distributions of larger number of particles, as we have described.

As before, we are going to find a weak solution to a system of the Laplace–Beltrami equations

$$\begin{aligned} \Delta _{\varOmega ^{(k)}_0} u_k = 0 \end{aligned}$$

defined on each particle domain $\varOmega ^{(k)}_0$, as in Fig. 5. On the internal boundaries of the domains, we set up continuity conditions, for example on $\partial \varGamma _3^{(2)}$ and $\partial \varGamma _2^{(3)}$

$$\begin{aligned} u_2\left. \right| _{\partial \varGamma _2^{(3)}}=u_3\left. \right| _{\partial \varGamma _3^{(2)}}=\mathrm{const}, \end{aligned}$$

(36)

$$\begin{aligned} \oint _{\partial \varGamma _2^{(3)}} \nabla u_2 \cdot \mathbf{n}_{s_2}|_{\partial \varGamma _{2}^{(3)}}\, \mathrm{d}l =-\oint _{\partial \varGamma _3^{(2)} }\nabla u_3 \cdot \mathbf{n}_{s_3}|_{\partial \varGamma _{3}^{(2)}}\, \mathrm{d}l. \end{aligned}$$

(37)

while on a few external boundaries, Dirichlet boundary conditions are set.

The numerical solution allows to calculate the total flux through the system by summing up the fluxes passing through the external contours, where the Dirichlet boundary conditions are set, either at the top of the pack or at the bottom using (10). The results are summarised in Fig. 8.

Remarkably, the reduced flux $Q_\mathrm{T}/ S_0 Q_0$ as a function of

$$\begin{aligned} \varPsi _0(H_B/2R_0)=\ln ^{-1} \left( \frac{1+ \sqrt{1-\left( \frac{H_B}{2R_0}\right) ^2} }{1- \sqrt{1-\left( \frac{H_B}{2R_0}\right) ^2}} \right) , \end{aligned}$$

(38)

where $H_B$ is the bridge size $H_B=2R_0\sin \theta _0$, behaves linearly in all configurations. This behaviour mirrors the flux dependence observed in azimuthally symmetric analytical solutions, see (15) or (31). The variations in the dependencies between different configurations are observed to be well within the statistical error expected in this case, error bar in Fig 8. At the same time, a comparison with a similar, but a regular arrangement, as in Fig. 3 at $R=R_0$ demonstrates that there is a clear cut contribution from the effects of tortuosity, solid line in Fig. 8.

Indeed, given identical porosity ($\phi \approx 50\%$) and mean particle size ($R/R_0=1$) in the regular, symmetric and randomly generated configurations, the normalised flux values differ by a factor of two $\tau _s\approx 2$, which can be regarded as the tortuosity of the surface flow. The value of the surface flow tortuosity $\tau _s\approx 2$ is consistent with the tortuosity values obtained in porous media in different conditions and configurations (Ghanbarian et al. 2013). For example, both hydraulic $\tau _h$ and diffusive $\tau _d$ tortuosities estimated in unsaturated porous media using different permeability models [often used in applications, for example Mualem (1976, 1978)] were found in between $1.5 \le \tau _{h,d}\le 2$ at $\phi =50\%$ (Ghanbarian et al. 2013).

It is important that the result, that is the ratio of the total flux in the random and regular configurations, does not practically depend on the size of the contour $H_B$, basically the size of the liquid bridge, and hence the value of saturation in the porous media. This implies that the observed effect is purely down to the distribution of contacts between the particles, but not the particular pathway on each single particle surface. That is, fundamentally, tortuosity in the surface diffusion processes is a geometric factor independent of the particular surface flow regime. At the same time, the pathways, on average, of course, does depend on the bridge size value $H_B$ leading to smaller permeability as the size of the contact area diminishes. This trend is expected, but essentially, the correction to the effective coefficient of diffusion

$$\begin{aligned} D(s)\propto \frac{K(s)}{(s-s_0)^{3/2}},\quad K(s)=\tau _s K_1(s) \end{aligned}$$

(39)

is only down to a single universal factor of two, $\tau _s$, representing the tortuosity effects in surface diffusion in particular porous media at low values of saturation, where $K_1(s)\propto |\ln (s-s_0)|^{-1}$ defined by (18) is the surface permeability of a single porous matrix element. Note, the value of $\tau _s\approx 2$ is also in agreement with experimental observations and a comparison of the superfast diffusion model with the data, where the tortuosity effects were estimated to reduce the effective permeability twofold (Lukyanov et al. 2019).

This is the main result of this study, which can be used in practical applications to calculate permeability in particular porous media. Basically, as the first step, one can calculate permeability of a single, representative element of the media or several elements to obtain some mean value and its dispersion. This way, permeability K(s), via (18), and the diffusion coefficient D(s) for the macroscopic model can be established in the first approximation. Macroscopic permeability K(s) or the diffusion coefficient D(s) then should be corrected by the universal factor of two, $\tau _s$, in the macroscopic diffusion model.

4 Conclusions

We have demonstrated that the Laplace–Beltrami method can be used to obtain permeability of particulate porous media at low saturation levels and to estimate contribution from the effects of tortuosity. Essentially, analytical results obtained using azimuthally oriented coupled particles can be used with a universal correcting prefactor $\tau _s$, as in (39), to estimate permeability of particle ensembles. That is, from the practical point of view, results obtained by analysing single representative element of particulate porous media can be translated into permeability of a particle composition.