## Abstract

A computational framework is developed for applying interfacial kinetic transport theory to predict water vapour permeability of porous media. Modified conservation equations furnish spatially periodic disturbances from which the average flux and, thus, the effective diffusivity is obtained. The equations are solved exactly for a model porous medium comprising parallel layers of gas and solid with arbitrary solid volume fraction. From the microscale effective diffusivity, a two-point boundary-value problem is solved at the macroscale to furnish the water vapour transport rate in membranes subjected to a finite RH differential. Then, the microscale model is implemented using a computational framework (extended finite-element method) to examine the role of particle size, aspect ratio and positioning for periodic arrays of aligned super-ellipses (model particles that pack with high density). We show that the transverse water vapour permeability can be reduced by an order of magnitude only when fibres with a high-aspect ratio cross section are packed in a periodic staggered configuration. Maximum permeability is achieved at intermediate micro-structural length scales, where gas-phase diffusion is enhanced by surface diffusion, but not limited by interfacial-exchange kinetics. The two-dimensional computations demonstrated here are intended to motivate further efforts to develop efficient computational solutions for realistic three-dimensional microstructures.

## 1. Introduction

Albaalbaki & Hill [1] recently set out a microscale model to quantify vapour permeability in porous media, specifically addressing the technologically important process of water vapour transmission in membranes and other microporous and nanoporous materials. This model was developed to predict water vapour transport in cellulose-fibre packing materials, but may be well suited to understanding moisture uptake and transport in composite materials and ceramics, particularly those with micro- and nano-dimensioned pores, which promote adsorption and surface diffusion fluxes. Other applications may be found in energy-storage and conversion devices. For example, Punckt *et al.* [2] highlight nanoscale porosity as enabling diffusive transport in batteries, fuel cells and solar cells. In cellulose-fibre-based barriers, the effective diffusivity depends on the adsorption of water molecules onto the solid surfaces and, possibly, absorption into the amorphous solid domains. The solid comprises closely packed nanofibrils (diameter approx. 10 nm) that are decorated with surface hydroxyl moieties, which promote multi-layer adsorption [3]. By solving coupled gas, solid and surface species conservation equations with interfacial boundary conditions that model molecular surface diffusion and interfacial exchange kinetics, the effective diffusivity may be predicted with knowledge of an equilibrium adsorption isotherm and solid-phase solubility; gas, solid and surface diffusion coefficients and interfacial kinetic rate constants.

Albaalbaki & Hill applied this model to a porous medium comprising randomly packed spheres, also adopting the self-consistent Maxwell approximation to predict the effective diffusivity at finite sphere solid volume fractions. In contrast to more conventional diffusion models, including the classical model of Maxwell [4], the extended Maxwell model [5] and parallel-pore diffusion models [6–8], the theory of Albaalbaki & Hill furnishes a permeability that can vary significantly with the microstructural length scale. Interestingly, because decreasing the pore size generally increases the specific interfacial surface area, one generally expects the flux from surface diffusion to increase with decreasing pore size. Therefore, decreasing the pore size is expected to increase the permeability. However, Albaalbaki & Hill found that the surface diffusion flux can become limited by interfacial exchange kinetics when the microstructural length scale is small. To help reduce water vapour transmission in hydrophilic materials, such as paper and ceramics, the theory suggests that the fibre or particle radius should be selected, when possible, according to the surface diffusivity, adsorption and absorption characteristics.

In this paper, we derive modified conservation equations and interfacial boundary conditions that can be implemented in a fully computational framework to rigorously address the role of microscale dimensionality and geometry. This establishes a well-defined boundary-value problem with periodic boundary conditions to compute effective diffusivities that can be applied in macroscale transport models of vapour transport in porous monoliths and membranes, for example. Thus, whereas Albaalbaki & Hill [1] solved their model for dilute sphere packings, the computational implementation undertaken here removes the geometrical approximations, enabling the impact of particle–particle and particle–matrix interactions to be fully explored. This general methodology has been adopted extensively in the literature to compute the hydrodynamic permeability or porous media [9]. In the hydrodynamic analogue, flow is driven by a mean pressure gradient or uniform body force, producing spatially periodic velocity and pressure disturbances. In stochastic microstructures, periodic artefacts from a finite computational domain vanish when the period is sufficiently large compared with the characteristic correlation length [10,11]. For diffusion in porous media, the analogue of a mean pressure gradient is the mean concentration gradient. The resulting periodic concentration perturbation to a (weak) average concentration gradient is the scalar product of the average gradient and a vector field that satisfies modified conservation equations and boundary conditions [12].

A microscale computational methodology is necessary to test analytical approximations, such as cell models and the Maxwell self-consistent approximations for random packings of spheres and discs. For dispersed porous media, computations explicitly account for disturbances to the average concentration gradient that arise from particle–particle interactions [13,14], which are absent in the theories based on a single-particle disturbance [4]. Note that cell models [12], which are appealing for the convenience with which they can furnish accurate predictions at low and high solid volume fractions [15,16], are also without rigorous justification. Finally, computations expedite studies exploring how the effective diffusivity varies between porous media with ordered, stochastic and anisotropic microstructures, which are beyond the capabilities of analytical calculation. Nowadays, it is possible to solve transport models in highly anisotropic and heterogeneous computational domains obtained from three-dimensional topography of real porous media [17,18], making it possible to predict permeability from empirical knowledge of the microstructural geometry.

To test the model, we solve the equations exactly for a lamellar porous medium and implement the microscale diffusivity at the macroscale to predict the water vapour transport rate (WVTR) in membranes with finite thickness when subjected to finite relative humidity (RH) differentials. While the lamellar microscale geometry is highly ideal, the calculations furnish valuable insights upon comparison to their two- and three-dimensional model counterparts: randomly packed spheres and discs. Note also that a lamellar microstructure furnishes an exact solution of the interfacial kinetic transport model for all volume fractions, which is valuable for testing computations, such as those undertaken in this paper using an extended finite-element methodology [19,20].

Accordingly, the model is applied at the microscale (representative elementary volume) for two-dimensional geometries. We report the effective diffusivity for ordered arrays of super-ellipses and discs, as prototypical geometries that achieve high solid volume fractions and anisotropy without particle overlap. For dispersed isotropic media, particle interactions play a weak role, as revealed by the self-consistent Maxwell approximation for randomly packed discs furnishing an accurate prediction for square arrays of discs at solid volume fractions up to the point of contact. For anisotropic microstructures, however, computations are necessary to reliably quantify the roles of solid volume fraction and particle aspect ratio, for example. We will draw comparisons to experiments in the literature, briefly exploring the roles of surface adsorption and surface mobility in determining overall permeability.

## 2. General microscale theory

Here, we formulate a general boundary value problem that can be solved with an arbitrary microscale geometry, as depicted schematically in figure 1*a*, to furnish the effective diffusivity in a periodic representative elementary control volume (shaded box in figure 1*a*) under conditions where a spatially uniform equilibrium state is perturbed from equilibrium by the application of a mean concentration gradient ** B**. By solving a sequence of such problems, each of which corresponds to a different local equilibrium, the flux on a macroscopic scale may be calculated with an effective diffusivity that depends on the local equilibrium state. This task is demonstrated in §4, after deriving an exact solution of the microscale problem for a lamellar porous medium in §3. A fully computational implementation of the microscale model is undertaken in §5 for rectangular and staggered arrays of high-aspect-ratio super-ellipses.

Albaalbaki & Hill [1] derived the following model equations that couple diffusion in solid and gas phases to interfacial kinetics and surface diffusion:
2.1
2.2
2.3
and
2.4Here, ** n** is an unit normal (from the solid to the gas phase); and

**∇**

_{i}=

**∇**−

**(**

*n***⋅**

*n***∇**) is the surface gradient operator. Equations (2.1) are steady diffusion equations for the gas and solid phases, equation (2.2) is a surface conservation equation, which balances surface diffusion with exchange fluxes between the interface and gas and solid phases, and equations (2.3) and (2.4) link the gas and solid diffusion fluxes to the interfacial exchange fluxes. The solid- and gas-phase concentration perturbations from the uniform equilibrium state,

*n*′

_{s}and

*n*′

_{g}, are scaled with the equilibrium gas-phase concentration

*n*

^{0}

_{g}, the surface concentration perturbation

*n*′

_{i}is scaled with the equilibrium surface concentration (

*K*is the partition coefficient from an adsorption isotherm), and lengths are scaled with a microstructural length scale ℓ, which Albaalbaki & Hill set equal to the sphere radius

*a*.

Note that measures the equilibrium solubility of the tracer in the solid phase; and *Π*-parameters are dimensionless model parameters involving diffusion coefficients (*D*_{s}, *D*_{g} and *D*_{i}) and interfacial kinetic exchange rate coefficients (*k*_{gi}, *k*_{gs}, *k*_{sg}, *k*_{gi} and *k*_{si}): , , ; and , , . Here, subscripts g, s and i denote gas, solid and interface, respectively. In this paper, the characteristic microstructural length scale for a lamellar microstructure is set to half the period of the minimum representative control volume, as shown in figure 1*b*; when we examine arrays of super-ellipses, we set ℓ to half the length of a square that has the same cross-sectional area.

We will seek numerical solutions of these equations when the concentration perturbations take (dimensional) forms *n*′_{g}=** B**⋅

**(1+**

*x**f*

_{g}),

*n*′

_{s}=

*H*

**⋅**

*B***(1+**

*x**f*

_{s}) and

*n*′

_{i}=

*K*

**⋅**

*B***(1+**

*x**f*

_{i}), where

*f*

_{g},

*f*

_{s}and

*f*

_{i}are dimensionless scalar functions of position in their respective domains

*Ω*

_{g},

*Ω*

_{s}and

*Γ*

_{i}, as depicted in figure 1

*a*; and

**is the (imposed) average gas-phase concentration gradient. The accompanying dimensionless perturbations are conveniently written [12,21] In this form, the perturbations are linear in**

*B***, and terms of the form**

*B***⋅**

*B***are periodic functions of position in a periodic representative control volume. This makes them very convenient to calculate in a computational environment, using discretization (e.g. finite-element, spectral or finite difference) or lattice methods (e.g. lattice-Boltzmann or moment propagation) for solving partial differential equations [22,23].**

*f*The Laplacians become ∇^{2}*n*′=∇^{2}** B**⋅

**, the normal gradients are**

*f***⋅**

*n***∇**(

**⋅**

*B***+**

*x***⋅**

*B***)=**

*f***⋅**

*n***+(**

*B***⋅**

*n***∇**)(

**⋅**

*B***), and the surface Laplacian is . Therefore, substituting these transformations into the foregoing conservation equations and interfacial boundary conditions furnishes 2.5 2.6 2.7 and 2.8where Note that these linear equations are independent of the magnitude of**

*f***, so one needs only to specify the orientation of**

*B***with respect to the microstructure to compute**

*B**g*

_{g}(

**),**

*x**g*

_{s}(

**) and**

*x**g*

_{i}(

**). From the solution of this boundary-value problem, the average flux (in dimensional variables) 2.9becomes where**

*x**Ω*=

*Ω*

_{g}+

*Ω*

_{s}(

**is the second-rank identity tensor), or 2.10where 2.11Note that the effective diffusivity tensor**

*δ*

*D*_{e}is defined by 2.12

## 3. One-dimensional solutions for a lamellar porous medium

Here, we consider the periodic unit cell in figure 1*b* with (dimensionless) coordinate |*x*|≤1. Note that the solid occupies |*x*|<*ϕ*, and the gas–solid interfaces are at |*x*|=*ϕ*. Linearity and symmetry considerations furnish an effective diffusivity tensor of the form
Thus, solving two relatively simple problems, each with ** B** directed either parallel or perpendicular to

**, furnish the scalar effective diffusivities**

*n**D*

^{⊥}

_{e}and . In the solid and gas phases, ∂

^{2}

*g*

_{g/s}/∂

*x*

^{2}=0, furnishing general solution

*g*

_{s}=

*B*(

*a*

_{s}

*y*+

*b*

_{s}) and when

**is perpendicular to**

*B***, and**

*n**g*

_{s}=

*B*(

*a*

_{s}

*x*+

*b*

_{s}) and when

**is parallel to**

*B***. Here, superscripts + and − identify the gas-phase solutions in the regions where**

*n**x*>

*ϕ*and

*x*<

*ϕ*, respectively.

### (a) Gradient perpendicular to the layers

When ** B** is perpendicular to the layers, ±

**⋅**

*B***=±**

*n**B*, (±

**⋅**

*n***∇**)=±∂⋅/∂

*x*and ∂⋅/∂

*y*=0, so the interfacial conservation equation requires with (since

*g*

^{±}

_{i}is periodic and finite). The boundary conditions (at

*x*=±

*ϕ*) require Note also that periodicity demands and . Therefore, there are seven equations with seven unknowns:

*a*

_{s},

*b*

_{s},

*a*

_{g},

*b*

^{±}

_{g},

*g*

^{±}

_{i}. Explicit solutions are provided in appendix A. Since the flux is uniform, the (dimensional) average flux can be written (in the gas phase) which is independent of the arbitrary value to which we assign

*b*

_{s}=0 (see appendix A). The effective diffusivity is therefore where an explicit formula for

*a*

_{g}is provided in appendix A. Note that the flux can also be written (in the solid phase) giving which is the same as calculated in the gas phase. There is no contribution of surface diffusion to the total flux, because there is no surface gradient in the direction of the bulk concentration gradient (

**⋅**

*τ***=0).**

*B*Figure 2 shows a representative solution of the foregoing boundary-value problem when ** B**⋅

**=**

*n**B*. Note that the perturbed concentration gradient in the gas phase has the opposite sign as the imposed average gradient

**. In the solid phase, however, the perturbation has the same sign as the average gradient. These gradients are analogous to the three-dimensional polarizations of a single solid sphere in an unbounded gas. As discussed in detail by Albaalbaki & Hill [1], the sign of this net polarization depends on the solubilities and phase diffusivities, ultimately determining whether the solid increases or decreases the effective diffusivity with respect to the gas-phase diffusivity.**

*B*### (b) Gradient parallel to the layers

When ** B** is directed parallel to the layers,

**⋅**

*B***=0, (**

*n***⋅**

*n***∇**)=0, ∂⋅/∂

*x*=0 and ∂⋅/∂

*y*=0 (since the solution must be periodic and there is no characteristic length scale in the

*y*-direction), so the interfacial conservation equation requires Therefore, (since

*g*

_{i}is periodic and finite); note also that , which demands (periodicity), so the boundary conditions (at

*x*=±

*ϕ*) require Again, setting

*b*

_{s}=0, the solution is

*g*

_{i}=

*b*

_{i}=0. Therefore, the (dimensional) average flux giving an effective diffusivity As expected, because there are no disturbances to the average gradients, the effective diffusivity is the volume-averaged diffusivity. Accordingly, the contribution of surface diffusion to the total flux is simply proportional to the specific interfacial surface area (

**⋅**

*τ***=ℓ**

*B*^{−1}

**), and the total flux is a monotonically decreasing function of the microstructural length scale ℓ.**

*B*The scaled effective diffusivities are plotted in figure 3 versus the solid volume fraction *ϕ* for several values of the solid-phase solubility parameter *H* with *Π*_{g}/*Π*_{s}=*D*_{g}/*D*_{s}=10^{4}. Here, the low solid-phase diffusivity relative to the gas-phase diffusivity furnishes *D*^{⊥}_{e}<*D*_{g} and *D*^{∥}_{e}<*D*_{g} unless *H* is large enough for the high solid-phase concentration to compensate for the low solid-phase mobility, i.e. .

## 4. Macroscale implementation

Following Albaalbaki & Hill [1], we calculated the macroscale flux through a membrane by numerically solving a two-point boundary problem where the local effective diffusivity varies according to the local gas-phase concentration . Here, *x* is the RH and *n**_{g} is the saturation concentration at the prevailing temperature and pressure. For these calculations, we adopted the well-known Guggenheim–Anderson–de Boer (GAB) theory adsorption isotherm [24] to specify
4.1which appears in the dimensionless parameter *Π*_{i}. Note that is the specific binding-site area, and the fractional occupancy
4.2Here, *c*≥0 is a dimensionless parameter that compares the water-binding site affinity with the water–water affinity, and *k* is a dimensionless parameter that allows the well-known Brunauer–Emmett–Teller (BET) theory [25] to better fit experimental data (by permitting a finite fractional occupancy when *x*=1). At room temperature, for example, Gupta & Chatterjee [8] measure *c*≈56 and *k*≈0.75 for paperboard (moisture content at monolayer coverage *m*=0.051 g water per gram solid), and Rhim & Lee [26] find *c*≈9.8 and *k*≈0.87 for Kraft paper (moisture content at monolayer coverage *m*=0.030 g water per gram solid). Note that the permeability *P* is defined by the macroscale flux ; it rigorously accounts for the tracer solubility and mobility through the interfacial kinetic model that furnishes the local effective diffusivity .

The permeability is often considered to be the product of a bulk solubility *S* and an effective (mobility) diffusivity *D*_{m}, i.e. *P*=*D*_{m}*S*. For example, in a porous medium where the concentration gradient is weak and, thus, *n*^{0}_{g}=const., *J*=*D*_{e}*B*=*D*_{m}*SB*, where *SB* is the gravimetric moisture content gradient. Letting *S*=*ϕH*+*sK*, where *s* is the surface area per unit volume, it follows that
4.3For arrays of non-overlapping discs, *s*=2*ϕ*/*a*; for arrays of non-overlapping spheres, *s*=3*ϕ*/*a* and for the lamellar microstructure above, *s*=ℓ^{−1}. Therefore, in a porous medium where the solid solubility is low (*ϕH*≪1) and the specific surface area is high, we find
when, for example, *ϕ*∼1, *θ*∼0.1, *a*_{m}∼0.1 Å , ℓ∼1 μm and .

Figure 4 shows the WVTR and dimensionless permeability *P*/*D*_{g} for a membrane with thickness *L*=100 μm when subjected to a total gas-phase concentration gradient *B*=*n**_{g}*Δ*RH/*L*, with *Δ*RH=0.5 at room temperature. In contrast to the calculations of Albaalbaki & Hill, this lamellar microstructure furnishes a permeability that decreases monotonically with increasing microstructural length scale. This is because the contribution of surface diffusion to the total flux increases with the specific interfacial surface area ℓ^{−1}. With the parameters adopted in this example, the limiting flux when becomes equal to the gas-phase flux, i.e. .

We calculated the gravimetric moisture content (GMC) as the mass of adsorbed and absorbed water per unit volume,
and report it as a mass percentage
where *ρ*=1400 kg m^{−3} is the (dry) solid-phase density, *M*_{w}≈0.018 kg mol^{−1} is the molecular weight of water and *N*_{A}≈6.022×10^{23} mol^{−1} is Avogardo’s number. Accordingly, figure 5 shows profiles of gas-phase moisture content (*a*) and gravimetric moisture content (*b*) for various relative humidities on each side of a membrane, for two values of the surface diffusivity. Despite different moisture contents, the WVTRs (reported in the figure caption) vary with the surface diffusivities, which vary here by four orders of magnitude; this reflects the small gas-phase volume fraction 1−*ϕ*=0.01 and the small microstructural length scale ℓ=100 nm, which together impart a large specific interfacial surface area. As discussed by Albaalbaki & Hill [1], WVTRs in cellulose fibre-based materials, e.g. paper, are generally because of the external gas-phase diffusion resistance, i.e. experiments typically measure the total resistance to moisture permeation, not the intrinsic resistance of the barrier membrane. Note also the absence of tortuously here: the lamellar pores transport vapour directly (uni-directionally) across the membrane.

## 5. Finite-element microscale implementation

To explore the roles of particle size, shape and packing configuration on permeability, we solved the model equations in two-dimensional doubly periodic arrays of super-ellipses. Such packings can achieve high solid volume fractions without particle overlap, and are therefore preferable to discs for mimicking dense wood-fibre packings, for example. Super-ellipses have a surface at positions *x* and *y* that satisfy
where *x*_{0} and *y*_{0} are the coordinates of the centre, the exponents *n*_{x} and *n*_{y} control the surface curvature, and *a*_{x} and *a*_{y} determine the particle size and aspect ratio. With a suitable choice of parameters, anisotropic porous media with ordered microstructure and high solid volume fraction can be constructed, without particle overlaps. All the computations in this paper were undertaken with *n*_{x}=6 and *n*_{y}=2, which produce the rounded ends and flat tops and bottoms, as seen in the figures below.

We validated our finite-element code, which we based on the extended finite-element methodology (XFEM) of Yvonnet *et al.* [19,20], by comparing the effective diffusivity of periodic square arrays or discs with an analytical solution of the model for random disc packings, according to the self-consistent Maxwell approximation, which we derived following Albaalbaki & Hill [1] (for sphere packings). Note that the XFEM is well suited to the tracer number density and flux being discontinuous across the gas-solid interface, since the elements containing the interface are enriched by extra degrees of freedom for , and . Results are presented in figure 6 at a solid volume fraction *ϕ*≈0.64 for a variety of the independent dimensionless parameters *Π*_{i} and *H*. The computations clearly capture the same qualitative changes in the effective diffusivity with respect to large changes in these parameters. Interestingly, there is excellent quantitative correspondence, even though there is no a priori reason why the analytical approximation and numerical solution should be so close at this moderately high solid volume fraction.

As a further test of the numerical fidelity, figure 7 compares effective diffusivities over a wide range of solid volume fractions. Here, the self-consistent Maxwell approximation for discs (solid line) is compared to finite-element computations for square arrays of discs and super-ellipses. As expected, the finite-element results for discs depart from the self-consistent approximation as the solid volume fraction approaches the close-packing limit *ϕ*≈0.78. On the other hand, because the super-ellipses do not overlap at higher solid volume fractions, the computations are close to the self-consistent approximation over a wide range of solid volume fractions. With unit aspect ratio, the particle shape and packing density evidently exert very weak influences on the permeability. For this reason, we now turn to highly anisotropic microstructures; namely rectangular and staggered arrays of super-ellipses.

For the computations reported below with two super-ellipses in the periodic unit cell, we specified the dimensionless width *L*_{x}ℓ^{−1} and height *L*_{y}ℓ^{−1} of the domain according to the particle aspect ratio *γ* and an area fraction parameter , which is the area fraction of rectangular inclusions having width 2*a*_{x} and height 2*a*_{y}. Accordingly,
It follows that
so changing the aspect ratio *γ* corresponds to area-conserving deformations of particles that have a characteristic size 2ℓ and cross-sectional area 4ℓ^{2}. Moreover, with specified and *γ*,

Colour maps of the concentration disturbances are shown in figures 8 and 9 with the mean concentration gradient ** B** directed in the [1,1] direction. Note that these are the linear superposition of the solutions with

**directed in the [1,0] and [0,1] directions, each of which specifies the effective diffusivities in the**

*B**x*- and

*y*-directions, respectively. Because we will examine the effective diffusivities when elongating the inclusions in the

*x*-direction (

*γ*>1), we refer to the accompanying effective diffusivities as and

*D*

^{⊥}

_{e}, respectively. The disturbances in figures 8 and 9 highlight how the particle aspect ratio and positioning can significantly affect the permeability. Here, the flux in the

*y*-direction is increasingly sensitive to the particle positioning as the aspect ratio increases. In the square configuration (panel

*a*), there are gas-phase channels in both directions, whereas the staggered configuration (panel

*b*) significantly increases the gas-phase tortuosity in the

*y*-direction, furnishing .

The influence of microstructural anisotropy on the permeability is quantified in figure 10 for rectangular (*a*) and staggered (** b**) configurations of high-aspect ratio particles at high packing density. In both panels, the permeability is plotted as a function of the microstructural length scale ℓ with fixed nominal volume fraction and an high aspect ratio

*γ*=10. For reference, the computational results (circles) are compared to three analytical theories:

*D*

^{∥}

_{e}for the lamellar microstructure (solid line);

*D*

_{e}for random sphere packing according to the self-consistent Maxwell approximation (blue dashed lines) [1]; and

*D*

_{e}for random disc packings (red dashed lines).

Recall, in lamellar porous media, *D*^{∥} (solid line) increases without bound as , because there are no gas-phase diffusion or interfacial exchange limitations to surface diffusion. In contrast to randomly packed spheres and discs, however, the high particle aspect ratio and alignment furnish an highly anisotropic permeability. In the staggered configuration (*b*), *D*^{⊥} is about an order of magnitude smaller than *D*^{⊥} over the entire range of particle sizes ℓ. This reflects the significantly higher gas-phase and surface tortuosity. In striking contrast, the rectangular particle placement (panel *a*) produces an isotropic permeability when transport is controlled by diffusion in the gas phase (large ℓ).

Similarly to the self-consistent Maxwell approximations (for random sphere and disc packings), the permeability is finite at asymptotically small and large ℓ. The various transitions in figure 10 are associated with ℓ controlling the magnitude and ratio of the dimensionless parameters: *Π*_{i}=*D*_{i}*K*/(ℓ^{2}*k*^{0}_{gi}), *Π*_{s}=*D*_{s}/(ℓ*k*^{0}_{gi}) and *Π*_{g}=*D*_{g}/(ℓ*k*^{0}_{gi}). Here, the adsorption isotherm and gas-phase kinetic-exchange coefficient *k*^{0}_{gi} furnish *Π*_{i}≈140ℓ^{−2}, *Π*_{g}≈100ℓ^{−1}, and *Π*_{s}≈0.001ℓ^{−1} with ℓ measured in nanometres. Thus, *Π*_{s}/*Π*_{g}=*D*_{s}/*D*_{g}=10^{−4} and *Π*_{i}/*Π*_{g}=*KD*_{i}/(*D*_{g}ℓ)≈1.4ℓ^{−1}. A physical interpretation of the results in figure 10 therefore rests on comparing gas-phase diffusion, surface diffusion and interfacial-kinetic mechanisms of transport.

We begin with the regime where . Here, and with . Thus, surface diffusion is weak compared to gas-phase diffusion, furnishing the large-ℓ plateaus in figure 10 where the flux is limited by diffusion in the continuous gas phase. Next, when ℓ∼100 nm, we find *Π*_{g}∼1 with and . In this regime, the extent to which surface diffusion can enhance the net flux depends on interfacial exchange kinetics. Similarly, when ℓ∼10 nm, we find *Π*_{i}∼1 and with *Π*_{i}/*Π*_{g}∼1. In this regime, interfacial kinetics limit the extent to which the flux can be enhanced by surface diffusion. Next, when ℓ∼3 nm, we have and with *Π*_{i}/*Π*_{g}∼1, showing that any enhancement of the net flux by surface diffusion is limited by interfacial exchange kinetics. Finally, when , we find and with . In this regime, the potential of surface diffusion to enhance the net flux is limited by interfacial exchange kinetics. Because the particles do not form a percolating network, the net flux again becomes limited by diffusion in the gas phase. Accordingly, the effective diffusivities plateau in the small-ℓ limit to the same value as in the large-ℓ limit.

Next, we briefly examine the concentration disturbances when the average gradient is directed perpendicular or parallel to the major particle axis. Figure 11 shows disturbances in the small-ℓ limit where, recall, surface diffusion is fast, but is limited by interfacial exchange kinetics and gas-phase diffusion. Panels 11*a*(i) and *b*(i) compare rectangular (*a*) and staggered (*b*) arrays with the mean concentration gradient directed along the *x*-axis. Note that the disturbances correspond to flux lines that transport vapour from the down-stream end of each fibre to the closest up-stream end; this occurs in a manner that minimizes the overall tortuosity, i.e. by maximizing the available surface diffusion conductance and minimizing the gas-phase diffusion resistance. In the staggered configuration, flux lines transport vapour across the mean gradient from the sides of each particle to the sides of its nearest neighbours. Panels 11*a*(ii) and *b*(ii) compare rectangular and staggered arrays with the mean concentration gradient directed along the *y*-axis. Again, the disturbances correspond to flux lines that transport vapour from one surface to another across the intervening gas-phase. While the concentration disturbances look qualitatively different, they both correspond to flux lines that minimize the overall tortuosity, subject to the low permeability of the solid, high conductivity of the surfaces and resistance of the intervening gas. Finally, panels 11*a*(iii) and *b*(iii) compare disturbances with the mean concentration gradient in the [1,1] direction. These are the linear superposition of the disturbances in the respective *a*(i),*b*(i) and *a*(ii),*b*(ii) panels.

The influence of particle aspect ratio on the longitudinal and perpendicular effective diffusivities, for super-ellipses in the staggered configuration, is examined in figure 12. Note that the longitudinal (perpendicular) diffusivity *D*^{∥} (*D*^{⊥}) increases (decreases) with *γ* because this decreases (increases) the tortuosity. Clearly, a high aspect ratio furnishes especially low effective diffusivities in the staggered configuration. This motivates aligning flattened fibres to minimize the vapour permeability.

The influence of RH is explored in figure 13. In figure 13*a*, *D*^{∥}/*D*_{g} and *D*^{⊥}/*D*_{g} for high-aspect-ratio super-ellipses (*γ*=10) in the staggered configuration are plotted versus the RH *x*. For convenient reference, the solid line is *D*^{∥}/*D*_{g} for a lamellar microstructure with the same ℓ=100 nm and . Note that with ℓ=100 nm, *γ*=10 and , it follows that and . Thus, the specific surface area available for surface diffusion is higher than in the lamellar porous medium with ℓ=100 nm. Note also that the actual solid volume fraction of the super-ellipse microstructure (*ϕ*≈0.82) is somewhat less than indicated by . Accordingly, the volume available for gas-phase diffusion is also higher than in the lamellar porous medium with *ϕ*=0.9. Based on these observations, the dashed line in each of figure 13*a*,*b* shows *D*^{∥}/*D*_{g} for a lamellar microstructure with ℓ=63 nm and *ϕ*=0.82. This furnishes a longitudinal diffusivity that is closer to the longitudinal diffusivity for super-ellipses, again highlighting the roles of surface- and gas-phase tortuosity.

In figure 13, we have adopted the GAB isotherm (*k*≠1) and a surface diffusivity *D*_{i} that varies with the fractional occupation according to the multi-layer adsorption theory of Chen & Yang [27]. Accordingly, the ratio of the surface diffusivity to the value at zero coverage is
5.1where 0≤*α*≤1 is the ratio of rate constants for the first and higher layers, and *θ*(*x*;*c*,*k*) is the fractional occupation (e.g. from the GAB or BET isotherm). We found that a diffusion coefficient that varies with the fractional occupancy was essential to furnish a dependence on the RH *x* that resembles literature data. This highlights the important role of surface diffusion in controlling moisture transport in porous media with fine microstructure.

Note that the effective diffusivity reflects the equilibrium partition coefficient , which enters the model via the dimensionless parameter *Π*_{i}=*KD*_{i}/(ℓ^{2}*k*^{0}_{gi}). With the GAB isotherm, the dimensionless partition coefficient is plotted in figure 13*c* versus the RH *x* for several values of *c* with *k*=0.95. While the fractional occupancy vanishes as , *K* has a minimum when the surface affinity is high (*c*>2); this minimum occurs at a RH that increases from *x*=0 to 0.5 as *c* increases from 2 to . When the surface diffusivity is a constant, *c* and *k* impact the effective diffusivity solely by the amount of vapour that adsorbs onto the solid and permeates the porous medium by surface diffusion. Similarly to the effective diffusivities , the mobility diffusion coefficients *D*^{∥/⊥}_{m} each have a maximum with respect to *x*. These reflect the minimum in *K* with respect to *x* (e.g. when *c*>2) as well as the maximum in the surface mobility, as measured by . As shown in figure 13*d*, the surface diffusivity *D*_{i} has a distinct maximum with respect to *x*. Note, however, that the general shape of these curves can change significantly according to the values of *c*, *k* and *α*. Here, *k*=0.95 and *α*=0.01 furnish maxima at relative humidities *x*>60% and the resulting changes in the effective diffusivities ( and *D*^{∥/⊥}_{m}) with respect to *x* bear close resemblance to literature measurements, as discussed below.

## 6. Discussion

From the temporal dynamics of the moisture uptake from paper samples subjected to step changes in the ambient RH, Bedane *et al.* [28] recently reported mobility diffusivities *D*_{m}(*x*), which they acquired from the local moisture content *m*=*Sn*^{0}_{g} satisfying a mass conservation equation ∂*m*/∂*t*=*D*_{m}∂^{2} m/∂*z*^{2}, where *D*_{m} can be assumed constant if the step change in the ambient RH is small compared with the initial RH (here, *z* is the position across the sample thickness, and *t* is the time). Similarly to the model predictions in figure 13*b*, their measured *D*_{m} (≈1×10^{−12}−8×10^{−12} m^{2} s^{−1}) also exhibit a maximum at relative humidities between 60 and 80%, with the lowest values occurring at the lowest *x*, and an intermediate value of *D*_{m} at the highest *x*. Note, however, that the general shape of the curves, similarly to our model, varies according to specific chemical and physical treatments (e.g. coating and calendaring), which change the fibre shape, volume fraction (*ϕ*) and surface properties (*c*, *k*, *α*).

In other experimental studies, effective diffusivities have been reported that increase monotonically, by a factor of about 5 as the RH increases from 0 to about 90% [8]; these changes are similar to the model predictions of *D*_{e}/*D*_{g} in figure 13*a*, provided we accept that *D*_{g}∼5×10^{−6} m^{2} s^{−1}, which is about five times smaller than the bulk diffusion coefficient of water in air (*D*_{g}≈2.5×10^{−5} m^{2} s^{−1} at *T*=25^{°}C).

Similarly, if we compare the model predictions of from figure 13*b* with *D*_{m}∼5×10^{−12} m^{2} s^{−1} from the experiments of Bedane *et al.* [28], then we find *D*_{g}∼3×10^{−6} m^{2} s^{−1}, which, again, is about one order of magnitude smaller than the continuum value for water vapour in air. The difference may be attributed to the shorter mean free path (i.e. Knudsen diffusion) in the nano-dimensioned spaces between wood fibres in paper. While the results in figure 13 have a nanometre-dimensioned fibre diameter ℓ∼100 nm, we found that increasing ℓ to 4 μm still furnishes . Note, however, that *D*^{⊥}_{e}/*D*_{g} is then ∼1, so if we take the experimental determination of *D*_{m}∼10^{−12} m^{2} s^{−1}, we again have (owing to Knudsen diffusion) giving *D*_{e}∼10^{−6} m^{2} s^{−1}.

Note that changing the fibre diameter demands a commensurate decrease in *a*_{m} to maintain the same equilibrium moisture content. Accordingly, increasing ℓ from 100 nm to 4 μm must be accompanied by a decrease in *a*_{m} from 1 to 0.025 Å . Such an unphysically high binding site density (*a*_{m} is less than the size of a water molecule) reinforces the notion that wood fibres have their own surface micro/nanostructure. This seems very plausible, since it is well known that the fibres are hierarchical assemblies of micro- and nanofibrillar cellulose.

This also poses some new and interesting questions for future research; specifically addressing how one should account for the widely disparate length scales in theoretical models. For example, it may be reasonable to first apply the model developed in this paper to the nanofibrillated surface microstructure where ℓ∼1 nm (nanofibril scale, rather than the wood fibre scale), and to then apply this effective diffusivity as the surface diffusivity for the macrofibre (ℓ∼1 μm), albeit with a renormalized interfacial kinetic transport model. This would require *Π*_{sg}=0 and a redefinition of *Π*_{i} in the Albaalbaki–Hill model; namely , where *D*′_{i} is now the effective diffusion coefficient for the porous surface layer, and , where is the equilibrium moisture content (same dimensions as *n*^{0}_{g}) in the surface layer, and *δ* is the layer thickness.

## 7. Concluding remarks

Our calculations show that low-vapour permeability, when there is significant surface adsorption and diffusion, can be achieved in (ideal) membranes comprising closely packed, staggered co-aligned fibres. While the permeability of these ideal, two-dimensional microstructures is sensitive to the microstructural length scale (e.g. fibre diameter), permeability is particularly sensitive to the volume fraction and geometrical tortuosity. This highlights the value of a computational methodology for complex microstructural geometries.

Nevertheless, the length-scale dependence of the permeability from finite-element computations is similar to the theoretical predictions of Albaalbaki & Hill [1] for the prototypical three-dimensional porous medium comprising randomly packed spheres. According to the self-consistent Maxwell approximation at large solid volume fractions, sphere (and disc) packings furnish moderate permeabilities at all length scales, with a maximum permeability occurring when interfacial kinetics limit surface diffusion.

Examining how the permeability depends on RH reveals how exquisitely sensitive the effective diffusivities *D*_{e} and *D*_{m} can be to the moisture adsorption isotherm (specifically the fractional occupancy) and surface diffusivity. Even with the few (*c*, *k* and *α*) parameters associated with the GAB isotherm and the multi-layer surface mobility theory of Chen & Yang [27], we found that a great variety of qualitatively different permeability dependencies on relative humidity can emerge. In this paper, we demonstrated qualitative correspondence with *c*=10, *k*=0.95 and *α*=0.01, the order of magnitudes of which are comparable to independent adsorption and surface diffusion studies (for cellulose fibre-based materials) in the literature.

We hope that the computational methodology set out in this paper will prompt further investigation, aided with the development of efficient computational algorithms to solve this model with more realistic three-dimensional microstructures. Although the finite-element methodology adopted in this study furnished an effective proof-of-concept (in two dimensions), significant improvements in computational performance may be required to tackle complex three-dimensional geometries in the manner that has been achieved in related studies of hydrodynamic permeability and conductivity using lattice-Boltzmann [17] and fast multi-pole methods [14], for example.

## Funding statement

Supported by the NSERC Innovative Green Wood Fibre Products Network.

## Appendix A. Coefficients for the lamellar porous medium

The following formulae were obtained by solving the linear algebraic equations in §3 of the main text using Mathematica. The solution is non-unique, depending on an arbitrary specification of *b*_{s}. Note, however, that this corresponds to the addition of an arbitrary constant to all the perturbations, which does not affect the resulting flux, as evident by *a*_{s} and *a*_{g} being independent of *b*_{s}.
A1
A2
A3
A4
A5
and
A6

- Received June 3, 2013.
- Accepted October 16, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.