Surface permeability of porous media particles and capillary transport

We have established previously, in a lead-in study, that the spreading of liquids in particulate porous media at low saturation levels, characteristically less than 10% of the void space, has very distinctive features in comparison to that at higher saturation levels. In particular, we have found that the dispersion process can be accurately described by a special class of partial differential equations, the super-fast non-linear diffusion equation. The results of mathematical modelling have demonstrated very good agreement with experimental observations. However, any enhancement of the accuracy and predictive power of the model, keeping in mind practical applications, requires the knowledge of the effective surface permeability of the constituent particles, which defines the global, macroscopic permeability of the particulate media. In the paper, we demonstrate how this quantity can be determined through the solution of the Laplace-Beltrami Dirichlet problem, we study this using the well-developed surface finite-element method.


Introduction
Liquid distributions and transport in particulate porous media, such as sand, at low saturation levels s, defined in our study as the ratio of the liquid volume V L to the volume of available voids V E in a sample volume element V , s = VL VE , have many distinctive features. Theoretically, as we have shown previously, the liquid dispersion can be described by a special class of mathematical models, the super-fast non-linear diffusion equation [1]. Unlike in the standard porous medium equation, which is a paradigm of research in porous media [2], in this special case, the non-linear coefficient of diffusion D(s) demonstrates divergent behaviour as a function of saturation s, D(s) ∝ (s − s 0 ) −3/2 , where s 0 is some minimal saturation level [1].
In practical applications, the analysis of this regime of wetting is crucial for studies of biological processes, such as microbial activity, and spreading of persistent (nonvolatile) liquids in soil compositions and dry porous media commonly found in arid natural environments and industrial installations [1,3].
If we consider liquid distributions on the grain size length scale, one would observe that when the saturation level s is reduced to (or below) the critical level s c ≈ 10%, the liquid domain predominantly consists of isolated liquid bridges formed at the point of particle contacts [1,[4][5][6][7][8], see fig. 1 for illustration. The formation of liquid bridges a e-mail: a.lukyanov@reading.ac.uk is characteristic for the so-called pendular regime of wetting [4][5][6][7][8][9]. In this regime, the liquid bridges are only connected via thin films formed on the rough particle surfaces and they serve as variable volume reservoirs, where the capillary pressure p depends directly on the amount of the liquid in the bridge V b Here, p 0 = 2γ R cos φ c , γ is the coefficient of the surface tension of the liquid, φ c is the contact angle made by the free surface of the liquid bridge with the rough solid surface of the constituent particles and R is an average radius of the porous medium particles [1,4,6]. The spreading process in such conditions only occurs over the rough surface of the elements of the particulate porous media connecting the liquid bridges, fig. 1.

Macroscopic formulation of the super-fast diffusion problem
Microscopically, the liquid creeping flow through the surface roughness of each particle can be described by a local Darcy-like relationship [10] between the surface flux density q and averaged (over some area containing many surface irregularities) pressure in the grooves ψ Here, μ is liquid viscosity and k m is the local coefficient of permeability of the rough surface, which proportional to the average amplitude of the surface roughness δ R , that is the width of the surface layer conducting the liquid flux, k m ∝ δ 2 R [10]. We note that, if the rough surface layer is not fully saturated with the liquid, parameter δ R should be interpreted as the characteristic width of the liquid layer within the rough surface layer. It is always assumed that δ R R, that is the amplitude of the surface roughness (or the width of the liquid layer) is always much smaller than the particle size.
Macroscopically, 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 a non-linear super-fast diffusion equation which directly follows from the conservation-of-mass principle ∂(φs) ∂t Here D 0 is the effective, macroscopic coefficient of nonlinear diffusion, s 0 is the minimal level of saturation, which can be only achieved when the liquid bridges cease to exist (s 0 ≈ 0.5%, see details in [1,7,8]), φ is porosity defined as φ = VE V , which is further assumed to be constant, and Q is the macroscopic flux density. The macroscopic flux density 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 Q · ndS, where n is the normal vector to the surface of the element.
Equation (3) can be obtained from (4) using (1), (2) and the spatial averaging theorem formulated in [11] assuming that [1]: -the rough surface area of the porous media particles is fully saturated with the liquid; -the liquid is incompressible; -the local Darcy's law (2) is observed on the rough particle surface elements. All three criteria are usually very well satisfied in practical applications, and we will further assume that this is the case. The approximation of the fully saturated rough surface layer is well fulfilled, if the characteristic pressure amplitude |ψ| is less than the capillary pressure amplitude defined on the length scale of the surface roughness δ R , which is of the order of δ R ∼ 1 μm in typical sands [12], as is demonstrated in [10]. That is, |ψ| < γ δR , and, for example for water (γ = 72 mN/m) at δ R = 1 μm, this results in |ψ| < 7.2 × 10 4 Pa. Otherwise, at larger absolute values of the (negative) capillary pressure, the liquid volume within the surface roughness layer would start to vary leading to variations of the effective liquid surface layer thickness δ R , though, it is not difficult to introduce a correction [3,[13][14][15]. Note, in the formulation (3), the effects of gravity were neglected assuming that the capillary length l c = γ/ρg 0 is much larger than the length scale associated with the gradient of the capillary pressure, that is l c √ δ R L 0 , where L 0 is the characteristic length scale of the wetting area. Here g 0 is the Earth gravity constant and ρ is the liquid density, so that for most liquids l c ∼ 1 mm. At the same time, taking δ R ≈ 1 μm and L 0 ≈ 10 mm, as it was in the experiments reported in [1], one gets  [1]. Here, parameter f φ = p0 2φ 3Nc 4π 1−φ φ , N c is a coordination number of the particles, that is the average number of contacts per a particle (in sands, typically, N c ≈ 7) and S e /S is the ratio of the effective area of entrances and exits of the liquid flow in a sample volume element with surface area S, see details in [1]. Note, that the ratio S e /S is defined in such a way, that the microscopic flux density q averaged over the liquid volume V l within a macroscopic sample volume element V , q l = V l q dV , if multiplied by the ratio q l Se S = Q, would result in the macroscopic average flux density Q.
The global surface permeability of the particles K is one of the main elements of the model that enables an accurate representation of the liquid dispersion at low saturation levels. On the other hand, this quantity is difficult to accurately estimate a priori. It is fully defined by the particle roughness and shape, and the dimension of the liquid bridge contact area, fig. 1. In this paper, we determine this important parameter on the basis of a solution to the Laplace-Beltrami problem in a representative case of a spherical (or nearly spherical) particle, which provides, as we will show, a reasonable approximation for the constituent elements of particulate porous media, such as sands.

Microscopic model of the surface permeability of the elements
Consider, as the simplest example, a spherical particle of radius R with a closed surface Γ , which is split into three sub-domains Ω 0 , Ω 1 and Ω 2 with the surface boundaries between them ∂Ω 1 and ∂Ω 2 , as is shown in fig. 2. The location of the sub-domains Ω 1 and Ω 2 to each other on the surface is fixed by the tilt angle α. The sub-domains Ω 1 and Ω 2 correspond to the contact area covered by the liquid in the bridges, while the surface flow, described by (2), takes place in Ω 0 .
Since the rough surface area of the particles is assumed to be fully saturated in creeping flow conditions [10], liquid pressure ψ, due to incompressibility of the liquid, should satisfy the Laplace-Beltrami equation defined on the surface of the sub-domain Ω 0 as follows from (2). Here, Δ Ω0 designates the Laplace-Beltrami operator, which is defined on the surface element Ω 0 through the surface gradient ∇ Ω0 tangential to the surface. Formally, let n Ω0 denote the unit normal to the surface Ω 0 then we define the surface gradient of ψ as ∇ Ω0 ψ := ∇ψ − (∇ψ · n Ω0 )n Ω0 and then the Laplace-Beltrami operator is defined as 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 δψ is negligible, that is δψ γ/δ R . 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 [9]. In the case when the surface layer is not fully saturated, parameter δ 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 to that in Ω 0 . So that, one can assume that which are the boundary conditions to the Laplace-Beltrami Dirichlet boundary value problem. The Dirichlet boundary value problem (5), (6) has a unique solution, which, if it is found, allows to calculate the total flux through the particle element where n is the normal vector to the domain boundaries ∂Ω 1,2 on the surface, δ 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 Ω 0 , for example the boundary ∂Ω 1 .
If the total flux Q 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 T ) if the flow is driven by the constant pressure difference ψ 2 − ψ 1 applied to the sides of the volume element.

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 ∂Ω 1 and ∂Ω 2 are oriented at the reflex angle α = π and have a circular shape. We use a spherical coordinate system with its origin at the particle centre and the polar angle θ counted from the axis of symmetry passing through the centre of the circular contour ∂Ω 1 . In this case, the Dirichlet boundary value problem (5), (6) admits an analytical solution, so that particle permeability can be determined explicitly. Indeed, problem (5), (6), if we assume that the liquid pressure ψ is a function of θ only and independent of the azimuthal angle, is equivalent to with the boundary conditions The analytic solution to problem (7), (8) after applying the boundary conditions can be represented in the following form: where Ψ 0 = 1 ln sin θ 1 sin θ 0 1 + cos θ 0 1 − cos θ 1 . One can now calculate the total flux So that, taking D = 2R, One can see that, if we take θ 1 = θ 0 , the permeability coefficient K 1 is divergent at θ 0 = π/2, as is expected, when the two contours move closer to each other and, at the same time, their radius R sin θ 0 increases, that is In the opposite limit, at θ 0 = 0, when the two contours move further away from each other and their radius decreases, the permeability coefficient tends to zero, that is Parametrically, the coefficient of permeability (10) 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 δ R , that is K 1 ∝ δ 3 R since it is anticipated that k m ∝ δ 2 R , so that evaluation of this parameter in applications is crucial for the accurate estimates of the liquid dispersion rates.
How does the result affect the super-fast diffusion model (3), 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 θ 1 = θ 0 , and, using an approximate relationship between the radius of curvature R sin θ 0 of the boundary contour ∂Ω 1 and the pendular ring volume [6], one can show As one can see from (11), the distinctive particle shape results in logarithmic correction to the main non-linear super-fast diffusion coefficient Apparently, the correction will mitigate to some extent the divergent nature of the dispersion at the very small saturation levels s ≈ s 0 , smoothing out the characteristic dispersion curves. Before we proceed to a general case, this would be instructive to consider, in qualitative terms, how specific is the permeability of spherical particles. We now compare coefficient of permeability (10) with the permeability of a cylinder of radius R sin θ 0 and length 2R with the same surface layer of thickness δ R . Such an element was often used in simple estimations of permeability in porous media [16]. It is not difficult to calculate the total flux through this element when there is a constant pressure difference (ψ 2 − ψ 1 ) applied to its ends where K c is the effective permeability of the cylindrical element.
One can observe, that in contrast to the case of spherical elements, the cylindrical approximation provides completely different correction to the non-linear coefficient of diffusion, if we presume similar scaling sin 2 θ 0 ≈ √ s − s 0 . Consider now a general case.

Surface permeability of a sphere in the case of arbitrary oriented boundaries
In the arbitrary case, when α = π, the Dirichlet boundary value problem (5), (6) does not possess known explicit solutions, so we make use of a classical surface finite-element technique introduced in [17]. See also [18] for an in depth review of state of the art innovations and uses pertaining to this class of method. Using this method we are able to numerically investigate the total flux and hence the permeability of the particle.
We begin by approximating the truncated surface element with a piecewise linear approximation through triangular elements, see fig. 3 for an example. In this setting, we are approximating the geometry with a polygon. This inherently introduces an error through the approximation Fig. 4. Verification of the numerical scheme in the azimuthally symmetric scenario. We plot the inverse mesh size against the error measured in the energy norm [17]. We observe the rate of convergence proven in [17] verifying that, asymptotically, the numerical approximation converges to the exact solution.
of the geometry. It is, however, well understood appearing as a "variational crime" [18]. We then discretise the Laplace-Beltrami operator over the polygon using piecewise linear finite elements. To test our numerical model we examine the azimuthally symmetric case, where the exact solution is known and given in (9). We then check convergence of the finite-element approximation to (9). The results are shown in fig. 4.
We make use of the numerical model generated to examine the dependency of the total flux, and hence the permeability of the truncated spherical element as a function of the tilt angle α, that is the position of the boundaries on the sphere at fixed values of the capillary pressure ψ 1 and ψ 2 . As in the azimuthally symmetric case, without much loss of generality, we consider circular boundaries. The size of the boundary contour, that is its radius R sin θ 0 (or R sin θ 1 ), will be characterized by the polar angle θ 0 (or θ 1 ) counted from the axis of symmetry of each contour and the particle radius R.

Results of numerical analysis and discussion
The distribution of pressure on the spherical surface is illustrated in fig. 5, while the typical total flux dependence on the tilt angle α is presented in fig. 6 at θ 0 = θ 1 and at fixed values of ψ 1 and ψ 2 . The distribution of pressure demonstrates relatively smooth variations in the range bounded by the prescribed boundary values, such that, as is expected in a diffusion problem, ψ 2 ≤ ψ ≤ ψ 1 . The value of the total liquid flux Q T through the spherical element decreases when the tilt angle increases and the boundary contours move further away from each other. At the same time, one readily observes, fig. 6, that at relatively large tilt angles, close to the reflex angle in the azimuthal symmetrical case, the total flux value and hence permeability of the surface elements, is close to that predicted on the basis of the azimuthally symmetric solution (10). This implies that the analytical result (10) and (11) can be  used in practical applications to obtain first-order corrections to the effective non-linear coefficient of dispersion in the super-fast diffusion model. One may notice that even at small tilt angles, when the two boundaries are located close to each other, one can still approximate the coefficient of permeability with the accuracy of 50%. We have verified numerically that in the general case the permeability coefficient of the particles demonstrates the same trends with variations of parameters θ 0 and θ 1 as in the azimuthally symmetric case.

Arbitrary particle shapes
Even low dispersed sand samples consist of grain particles, which are only approximately spherical [12]. Therefore, we consider arbitrary surface elements obtained by perturbations of a sphere preserving surface smoothness. Based on our methodology, we examine numerical solutions to the Laplace-Beltrami Dirichlet boundary value problem (5) set on such perturbed particle surfaces to calculate the total volumetric flux, which is the measure of the surface permeability. To separate the effects of the particle shape from the effects of the boundary shape on the particle surface permeability and for the sake of comparison with the permeability of spherical particles, we consider circular boundary contours oriented to each other as in the azimuthally symmetric case, fig. 7. The size of the boundary contour, that is its radius R sin θ 0 (or R sin θ 1 ), will be characterized by the polar angle θ 0 (or θ 1 ) counted from the axis of symmetry of each contour and the radius of the sphere used to obtain the perturbed surface element R. The first particle shape, we have examined, is shown in fig. 7 with the distribution of the liquid pressure indicated by the colour map. For the sake of comparison, we have chosen the same boundary conditions as in in the case of spherical shapes, that is ψ 1 /ψ c = 0.8 and ψ 2 /ψ c = 0.2, with the same contour sizes, that is θ 1 = θ 0 = π/8 oriented at α = π. As is expected, the total volumetric flux, in this case Q T , is reduced in comparison with that, Q 0 , through the spherical particle shape Q T /Q 0 ≈ 0.86, since some pathways connecting two boundary contours became much longer, as one can see from fig. 7. Despite, at first glance, strong variations of the original spherical shape, the observed effect is not dramatic and is on the scale of the change of the surface area demonstrating that the spherical shape provides a good approximation in general to obtain estimates of the surface permeability. Indeed, the total increase of the surface area due to the perturbation was S a /S 0 ≈ 1.2, where S 0 = 4πR 2 cos θ 0 is the surface area of the original truncated spherical particle, so that the characteristic size of the particle calculated via R a = R S a /S 0 ≈ 1.1 R. We note though that the actual parameter defining the particle permeability is expected to be an effective length of the pathways connecting the boundary contours.
In general, the effective pathway length scale is not so easy to estimate, therefore, to understand the role of this effective parameter, consider now specific systematic changes of the original spherical shape of radius R via a transformation of the form where θ and φ are the polar and azimuthal angles of the spherical coordinate system. The obtained surface profile is demonstrated in fig. 8 at A s = 0.15 and m = n = 5. As in the previous case, the boundary contours are circular, identical (θ 1 = θ 0 ) and are not perturbed. The smoothness of the perturbed surface shape was achieved via a spline approximation at the boundary contours during the mesh generation and further refinement of the mesh. In this procedure, a smooth surface profile is created with two small boundary regions, which are not exactly described by the transformation (12). In what follows, we fix the parameters of the perturbation transformation m = n = 5 and consider only variations of the amplitude A s . Variation of the total flux Q T through such elements with the amplitude of the perturbation A s is shown in fig. 9.
The characteristic arc length L p of the perturbed shape can be estimated by means of Fig. 9. Non-dimensional total flux QT /Q0 as a function of the shape perturbation amplitude A 2 s at fixed values of ψ1/ψc = 0.8 and ψ 2/ψc = 0.2, α = π, m = n = 5 and θ0 = θ1 = π/8. Here Q 0 is the total flux value through the unperturbed spherical element in similar conditions. The numerical results obtained at medium resolution (maximum mesh size h/R ≈ 0.03) are shown by symbols and the solid line is the best fit Q T /Q0 = 1 − CsA 2 s at Cs = 1.7. The approximation error is about the symbol size.
at A 2 s m 2 /8 1 and θ 0 = θ 1 . The estimate follows from the definition of L p along the meridian line (φ = const) taking into account that A 2 s m 2 /8 1 and applying averaging in the azimuthal direction, that is over φ, That is after averaging over the azimuthal angle and neglecting the contribution of the term of the order of sin θ 0 /2m 1 Since the total volumetric flux is expected to be proportional to the pressure gradient, one can anticipate that its dependence on the effective arc length would follow Q T /Q 0 ≈ (π−2θ0)R Lp ≈ 1 − A 2 s m 2 8 . As one can see, fig. 9, the numerically calculated total flux dependence does follow the trend suggested by scaling of the arc length L p , the match though is not perfect. We found from the best fit Q T /Q 0 = 1 − C s A 2 s , fig. 9, that C s = 1.7, while the value C s ≈ 3.1 would be expected. This implies that the surface diffusion over uneven landscapes is a slightly more complex phenomenon than that one would expect from the simple scaling suggested by the effective pathways length. We note, in that respect, that the methodology and the numerical treatment of the Laplace-Beltrami problem developed are particularly indispensable, where there is no simple way of estimating the effective parameter L p , for example over strongly heterogeneous surface profiles with large areas inaccessible to the liquid flow.